0001 function fwd_model= show_pseudosection( fwd_model, data, orientation)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019 if ischar(fwd_model) && strcmp(fwd_model,'UNIT_TEST'); do_unit_test; return; end
0020
0021 if ~isfield(fwd_model, 'misc') || ~isfield(fwd_model.misc,'sizepoint')
0022 if size(fwd_model.electrode,2) <= 16
0023 fwd_model.misc.sizepoint= 500;
0024 elseif size(fwd_model.electrode,2) <= 32
0025 fwd_model.misc.sizepoint= 200;
0026 else
0027 fwd_model.misc.sizepoint= 50;
0028 end
0029 end
0030
0031 if nargin<3
0032 try
0033 orientation= fwd_model.show_pseudosection.orientation;
0034 catch
0035 orientation = 'horizontaldownward';
0036 end
0037 end
0038
0039 if ~isfield(fwd_model,'show_pseudosection') || ~isfield(fwd_model.show_pseudosection,'fs')
0040 fwd_model.show_pseudosection.fs=16;
0041 end
0042
0043 if iscell(orientation)
0044 if strcmp(orientation{2},'yz')
0045 fwd_model.nodes= fwd_model.nodes(:,[2 3 1]);
0046 rotation= 'yz';
0047 elseif strcmp(orientation{2},'xz')
0048 fwd_model.nodes= fwd_model.nodes(:,[1 3 2]);
0049 rotation= 'xz';
0050 end
0051 orientation= orientation{1};
0052 end
0053 location34= [];
0054 electrodesLocation= [];
0055 switch(upper(orientation))
0056 case 'HORIZONTALDOWNWARD'; [depth,location]= plotPseudoSectionProfileDown(fwd_model,data);
0057 case 'VERTICAL'; [depth,location]= plotPseudoSectionProfileVert(fwd_model,data);
0058 case 'SPECIFICHORIZONTALDOWNWARD'; [depth,location,location34,electrodesLocation]= plotPseudoSectionProfileSpecificDown(fwd_model,data);
0059 case 'SPECIFICHORIZONTALUPWARD'; [depth,location,L,d]= plotPseudoSectionProfileSpecificUp(fwd_model,data);
0060 case 'CIRCULAROUTSIDE'; [depth,location]= plotPseudoSectionCircularOut(fwd_model,data);
0061 case 'CIRCULARINSIDE'; [depth,location]= plotPseudoSectionCircularIn(fwd_model,data);
0062 case 'CIRCULARVOLCANO'; [depth,location]= plotPseudoSectionCircularInVolcano(fwd_model,data);
0063 case 'HORIZONTALLAYERS'; [depth,location,electrodesLocation]= plotPseudoSectionHorizontalLayers(fwd_model,data);
0064 otherwise;
0065 error('No orientation of type "%s" available', upper(orientation));
0066 end
0067
0068 if exist('rotation','var')
0069 if strcmp(rotation,'yz')
0070 fwd_model.nodes= fwd_model.nodes(:,[3 1 2]);
0071 elseif strcmp(rotation,'xz')
0072 fwd_model.nodes= fwd_model.nodes(:,[1 2 3]);
0073 end
0074 end
0075
0076 fwd_model.show_pseudosection.depth= depth;
0077 fwd_model.show_pseudosection.location= location;
0078 fwd_model.show_pseudosection.location34= location34;
0079 fwd_model.show_pseudosection.electrodesLocation= electrodesLocation;
0080
0081
0082 end
0083
0084 function [depth,location,electrodesLocation]= plotPseudoSectionHorizontalLayers(fmdl,data)
0085 fs= fmdl.show_pseudosection.fs;
0086 depthPrecision= fmdl.show_pseudosection.depthPrecision;
0087 depthRatio= fmdl.show_pseudosection.depthRatio;
0088
0089 xprecision= fmdl.show_pseudosection.xprecision;
0090 yprecision= fmdl.show_pseudosection.yprecision;
0091 locFun= fmdl.show_pseudosection.locFun;
0092 locVar= fmdl.show_pseudosection.locVar;
0093
0094 depthFun= fmdl.show_pseudosection.depthFun;
0095 depthVar= fmdl.show_pseudosection.depthVar;
0096
0097 if ~isfield(fmdl.show_pseudosection,'uniqueDir');
0098 uniqueDir= 'x';
0099 else
0100 uniqueDir= fmdl.show_pseudosection.uniqueDir;
0101 end
0102 [elec_posn,elecNumber]= electrodesPosition(fmdl);
0103 if isfield(fmdl,'misc') && isfield(fmdl.misc,'elec_posn')
0104 elec_posn= fmdl.misc.elec_posn;
0105 end
0106 me= (mod(elec_posn,1));
0107 me(me==0)= 1;
0108 precision_position= min(me(:));
0109 xposition_elec= reshape(elec_posn(elecNumber,1),[],4);
0110 yposition_elec= reshape(elec_posn(elecNumber,2),[],4);
0111 zposition_elec= reshape(elec_posn(elecNumber,3),[],4);
0112 zposition_elec= zposition_elec-min(zposition_elec(:));
0113
0114 electrodesLocation= xposition_elec;
0115 if isfield(fmdl.show_pseudosection,'elecsUsed')
0116 elecsUsed= fmdl.show_pseudosection.elecsUsed;
0117 else
0118 elecsUsed= 1:4;
0119 end
0120
0121 location12x= mean(xposition_elec(:,elecsUsed(1:2)),2);
0122 location34x= mean(xposition_elec(:,elecsUsed(3:4)),2);
0123 location1234x= mean(xposition_elec(:,elecsUsed(1:4)),2);
0124 location12y= mean(yposition_elec(:,elecsUsed(1:2)),2);
0125 location34y= mean(yposition_elec(:,elecsUsed(3:4)),2);
0126 location1234y= mean(yposition_elec(:,elecsUsed(1:4)),2);
0127
0128 depth12x= diff(xposition_elec(:,elecsUsed(1:2)),1,2);
0129 depth34x= diff(xposition_elec(:,elecsUsed(3:4)),1,2);
0130 depth12y= diff(yposition_elec(:,elecsUsed(1:2)),1,2);
0131 depth34y= diff(yposition_elec(:,elecsUsed(3:4)),1,2);
0132
0133 x= cell(length(locVar),1);
0134 for i= 1:length(locVar)
0135 x{i}= eval(locVar{i});
0136 end
0137 location= locFun(x);
0138
0139 y= cell(length(depthVar),1);
0140 for i= 1:length(depthVar)
0141 y{i}= eval(depthVar{i});
0142 end
0143 depth= depthFun(y);
0144 [LU,is,js]= unique(location);
0145 dL= min(diff(LU));
0146 for j= 1:length(is)
0147 locationis= location(js==j);
0148 depthis= depth(js==j);
0149 if length(unique(depthis)) ~=length(depthis)
0150 [depthisU,isl,jsl]= unique(depthis);
0151 for k= 1:length(isl)
0152 if length(find(jsl==k))>1
0153 shift= (linspace(-dL/2,dL/2,length(find(jsl==k))))';
0154 locationis(jsl==k)= locationis(jsl==k)+shift;
0155 end
0156 end
0157 location(js==j)= locationis;
0158 end
0159 end
0160
0161 P= [depth location];
0162 Pu= unique(P,'rows','stable');
0163 if size(P,1)~= size(Pu,1)
0164 disp(' ')
0165 disp('P not unique !!!')
0166 disp(' ')
0167 keyboard
0168 end
0169
0170 scatter(location,depth,fmdl.misc.sizepoint,data,'filled','MarkerEdgeColor','k');
0171 hold on;
0172 xlabel('Distance (m)','fontsize',fs,'fontname','Times','interpreter','latex');
0173 ylabel('Distance (m)','fontsize',fs,'fontname','Times','interpreter','latex')
0174 axis equal; axis tight;
0175 bx= get(gca,'xlim'); by= get(gca,'ylim');
0176 set(gca,'fontsize',fs,'fontname','Times','dataAspectRatio',[1 (by(2)-by(1))/(bx(2)-bx(1))*3 1])
0177 end
0178
0179
0180 function [depth,location,location34,electrodesLocation]= plotPseudoSectionProfileSpecificDown(fmdl,data)
0181 fs= fmdl.show_pseudosection.fs;
0182 [elec_posn,elecNumber]= electrodesPosition(fmdl);
0183 xposition_elec= reshape(elec_posn(elecNumber,1),[],4);
0184 yposition_elec= reshape(elec_posn(elecNumber,2),[],4);
0185 zposition_elec= reshape(elec_posn(elecNumber,3),[],4);
0186 zposition_elec= zposition_elec-min(zposition_elec(:));
0187 if isfield(fmdl.show_pseudosection,'minx') && isfield(fmdl.show_pseudosection,'maxx') && ...
0188 isfield(fmdl.show_pseudosection,'miny') && isfield(fmdl.show_pseudosection,'maxy')
0189 minx= fmdl.show_pseudosection.minx;
0190 maxx= fmdl.show_pseudosection.maxx;
0191 miny= fmdl.show_pseudosection.miny;
0192 maxy= fmdl.show_pseudosection.maxy;
0193 else
0194 xposition_elec= xposition_elec-min(xposition_elec(:));
0195 yposition_elec= yposition_elec-min(yposition_elec(:));
0196 minx= min(xposition_elec(:));
0197 maxx= max(xposition_elec(:));
0198 miny= min(yposition_elec(:));
0199 maxy= max(yposition_elec(:));
0200 end
0201
0202 if isfield(fmdl.show_pseudosection,'xToDeduce') && isfield(fmdl.show_pseudosection,'yToDeduce')
0203 if strcmp(fmdl.show_pseudosection.xToDeduce,'minx')
0204 xToDeduce= minx;
0205 elseif strcmp(fmdl.show_pseudosection.xToDeduce,'maxx')
0206 xToDeduce= maxx;
0207 end
0208 if strcmp(fmdl.show_pseudosection.yToDeduce,'miny')
0209 yToDeduce= miny;
0210 elseif strcmp(fmdl.show_pseudosection.yToDeduce,'maxy')
0211 yToDeduce= maxy;
0212 end
0213 else
0214 xToDeduce= minx;
0215 yToDeduce= miny;
0216 end
0217 rposition_elec= sqrt((xposition_elec-xToDeduce).^2 + ...
0218 (yposition_elec-yToDeduce).^2);
0219
0220 if isfield(fmdl.show_pseudosection,'elecsUsed')
0221 elecsUsed= fmdl.show_pseudosection.elecsUsed;
0222 else
0223 elecsUsed= 3:4;
0224 end
0225
0226 if isfield(fmdl.show_pseudosection,'dirAxis')
0227 dirAxis= fmdl.show_pseudosection.dirAxis;
0228 else
0229 dirAxis= 'X';
0230 end
0231
0232 if isfield(fmdl.show_pseudosection,'depthRatio')
0233 depthRatio= fmdl.show_pseudosection.depthRatio;
0234 else
0235 depthRatio= 3;
0236 end
0237
0238 if isfield(fmdl.show_pseudosection,'depthPrecision')
0239 depthPrecision= fmdl.show_pseudosection.depthPrecision;
0240 else
0241 depthPrecision= 1;
0242 end
0243
0244 if isfield(fmdl.show_pseudosection,'depthLevel')
0245 depthLevel= fmdl.show_pseudosection.depthLevel;
0246 else
0247 depthLevel= 0;
0248 end
0249
0250 switch dirAxis
0251 case 'X'
0252 location= mean(xposition_elec(:,elecsUsed),2);
0253 depth12= -abs(xposition_elec(:,elecsUsed(1))-xposition_elec(:,elecsUsed(2)))/depthRatio;
0254 depth34= -abs(xposition_elec(:,elecsUsed(3))-xposition_elec(:,elecsUsed(4)))/depthRatio;
0255 electrodesLocation= unique(round(xposition_elec/depthPrecision)*depthPrecision);
0256 case 'Y'
0257 location= mean(yposition_elec(:,elecsUsed),2);
0258 depth= -abs(yposition_elec(:,elecsUsed(1))-yposition_elec(:,elecsUsed(2)))/depthRatio;
0259 electrodesLocation= unique(round(yposition_elec/depthPrecision)*depthPrecision);
0260 case 'Z'
0261 location= mean(zposition_elec(:,elecsUsed),2);
0262 depth= -abs(zposition_elec(:,elecsUsed(1))-zposition_elec(:,elecsUsed(2)))/depthRatio;
0263 electrodesLocation= unique(round(zposition_elec/depthPrecision)*depthPrecision);
0264 case 'R'
0265 location= mean(rposition_elec(:,elecsUsed),2);
0266 location12= mean(rposition_elec(:,elecsUsed(1:2)),2);
0267 location34= mean(rposition_elec(:,elecsUsed(3:4)),2);
0268 depth12= abs(rposition_elec(:,elecsUsed(1))-rposition_elec(:,elecsUsed(2)));
0269 depth34= abs(rposition_elec(:,elecsUsed(3))-rposition_elec(:,elecsUsed(4)));
0270 electrodesLocation= rposition_elec;
0271 end
0272
0273 precision= floor(min(diff(unique(rposition_elec(:)))));
0274 if precision==0
0275 precision= ceil(min(diff(unique(rposition_elec(:)))));
0276 end
0277
0278 depth12= round(depth12/precision)*precision;
0279 depth34= round(depth34/precision)*precision;
0280 location12= round(location12/precision)*precision;
0281 location34= round(location34/precision)*precision;
0282
0283 miSdepth12= min(diff(unique(depth12)));
0284 if ~isempty(miSdepth12)
0285 depth34reS= (depth34-min(depth34))/max(depth34)*miSdepth12;
0286 else
0287 depth34reS= 0;
0288 end
0289 depth= depth12 +depth34reS;
0290 depth= round(-(depth)/depthRatio/precision)*precision;
0291
0292 [depthU,is,js]= unique(depth);
0293 for j= 1:length(is)
0294 changeDepth= 0;
0295 clear locationis depthis
0296 location34is= location34(js==j);
0297 location12is= location12(js==j);
0298 depthis= depth(js==j);
0299 datais= data(js==j);
0300 if length(unique(location34is))==length(location34is)
0301 locationis= location34is;
0302 else
0303 [location34isU,isl,jsl]= unique(location34is);
0304 for k= 1:length(isl)
0305 shift= (location12is(jsl==k)-mean(location12is(jsl==k)))/(min(depth34));
0306 locationis(jsl==k)= round((location34is(jsl==k)+shift)/precision)*precision;
0307 end
0308 if length(unique(locationis))~=length(locationis)
0309 depthisShifted= depth(js==j)*0;
0310 for k= 1:length(isl)
0311 shift= ((location12is(jsl==k)-mean(location12is(jsl==k)))/(min(depth34)*2));
0312 if length(shift)==2 && round(shift(1)/precision)*precision==0
0313 shift= ((location12is(jsl==k)-mean(location12is(jsl==k)))/(min(depth34)*2));
0314 end
0315 if length(shift)==3
0316 xtremeShift= min([abs(min(shift)) abs(max(shift))]);
0317 if round(xtremeShift(1)/precision)*precision==0
0318 xtremeShift= max([abs(min(shift)) abs(max(shift))]);
0319 end
0320 shift= [-xtremeShift 0 xtremeShift]';
0321 end
0322 depthisShifted(jsl==k)= depthis(jsl==k)+shift;
0323 changeDepth= 1;
0324 end
0325 end
0326 end
0327 location(js==j)= locationis;
0328 if changeDepth
0329 depth(js==j)= depthisShifted;
0330 end
0331 end
0332
0333 depth= depth+depthLevel;
0334 scatter(location,depth,fmdl.misc.sizepoint,data,'filled','MarkerEdgeColor','k');
0335 hold on;
0336 xlabel('Distance (m)','fontsize',fs,'fontname','Times','interpreter','latex');
0337 ylabel('Pseudo-depth (m)','fontsize',fs,'fontname','Times','interpreter','latex')
0338 axis equal; axis tight;
0339 bx= get(gca,'xlim'); by= get(gca,'ylim');
0340 set(gca,'fontsize',fs,'fontname','Times','dataAspectRatio',[1 (by(2)-by(1))/(bx(2)-bx(1))*3 1])
0341 end
0342
0343 function [depth,location,L,d]= plotPseudoSectionProfileSpecificUp(fmdl,data)
0344 fs= fmdl.show_pseudosection.fs;
0345 [elec_posn,elecNumber]= electrodesPosition(fmdl);
0346 xposition_elec= reshape(elec_posn(elecNumber,1),[],4);
0347 yposition_elec= reshape(elec_posn(elecNumber,2),[],4);
0348 zposition_elec= reshape(elec_posn(elecNumber,3),[],4); zposition_elec= zposition_elec-min(zposition_elec(:));
0349 if isfield(fmdl.show_pseudosection,'minx') && isfield(fmdl.show_pseudosection,'maxx') && ...
0350 isfield(fmdl.show_pseudosection,'miny') && isfield(fmdl.show_pseudosection,'maxy')
0351 minx= fmdl.show_pseudosection.minx;
0352 maxx= fmdl.show_pseudosection.maxx;
0353 miny= fmdl.show_pseudosection.miny;
0354 maxy= fmdl.show_pseudosection.maxy;
0355 else
0356 xposition_elec= xposition_elec-min(xposition_elec(:));
0357 yposition_elec= yposition_elec-min(yposition_elec(:));
0358 minx= min(xposition_elec(:));
0359 maxx= max(xposition_elec(:));
0360 miny= min(yposition_elec(:));
0361 maxy= max(yposition_elec(:));
0362 end
0363
0364 if isfield(fmdl.show_pseudosection,'xToDeduce') && isfield(fmdl.show_pseudosection,'yToDeduce')
0365 if strcmp(fmdl.show_pseudosection.xToDeduce,'minx')
0366 xToDeduce= minx;
0367 elseif strcmp(fmdl.show_pseudosection.xToDeduce,'maxx')
0368 xToDeduce= maxx;
0369 end
0370 if strcmp(fmdl.show_pseudosection.yToDeduce,'miny')
0371 yToDeduce= miny;
0372 elseif strcmp(fmdl.show_pseudosection.yToDeduce,'maxy')
0373 yToDeduce= maxy;
0374 end
0375 else
0376 xToDeduce= minx;
0377 yToDeduce= miny;
0378 end
0379 rposition_elec= sqrt((xposition_elec-xToDeduce).^2 + ...
0380 (yposition_elec-yToDeduce).^2);
0381
0382
0383
0384
0385 if isfield(fmdl.show_pseudosection,'elecsUsed')
0386 elecsUsed= fmdl.show_pseudosection.elecsUsed;
0387 else
0388 elecsUsed= 3:4;
0389 end
0390
0391 if isfield(fmdl.show_pseudosection,'dirAxis')
0392 dirAxis= fmdl.show_pseudosection.dirAxis;
0393 else
0394 dirAxis= 'X';
0395 end
0396
0397 if isfield(fmdl.show_pseudosection,'depthRatio')
0398 depthRatio= fmdl.show_pseudosection.depthRatio;
0399 else
0400 depthRatio= 3;
0401 end
0402
0403 if isfield(fmdl.show_pseudosection,'depthPrecision')
0404 depthPrecision= fmdl.show_pseudosection.depthPrecision;
0405 else
0406 depthPrecision= 1;
0407 end
0408
0409 if isfield(fmdl.show_pseudosection,'depthLevel')
0410 depthLevel= fmdl.show_pseudosection.depthLevel;
0411 else
0412 depthLevel= 0;
0413 end
0414
0415 switch dirAxis
0416 case 'X'
0417 location= mean(xposition_elec(:,elecsUsed),2);
0418 depth= abs(xposition_elec(:,elecsUsed(1))-xposition_elec(:,elecsUsed(2)))/depthRatio;
0419 electrodesLocation= unique(round(xposition_elec/depthPrecision)*depthPrecision);
0420 case 'Y'
0421 location= mean(yposition_elec(:,elecsUsed),2);
0422 depth= abs(yposition_elec(:,elecsUsed(1))-yposition_elec(:,elecsUsed(2)))/depthRatio;
0423 electrodesLocation= unique(round(yposition_elec/depthPrecision)*depthPrecision);
0424 case 'Z'
0425 location= mean(zposition_elec(:,elecsUsed),2);
0426 depth= abs(zposition_elec(:,elecsUsed(1))-zposition_elec(:,elecsUsed(2)))/depthRatio;
0427 electrodesLocation= unique(round(zposition_elec/depthPrecision)*depthPrecision);
0428 case 'R'
0429 location= mean(rposition_elec(:,elecsUsed),2);
0430 depth= abs(rposition_elec(:,elecsUsed(1))-rposition_elec(:,elecsUsed(2)))/depthRatio;
0431 electrodesLocation= unique(round(rposition_elec/depthPrecision)*depthPrecision);
0432 end
0433 L= (rposition_elec(:,2)-rposition_elec(:,1))/2;
0434 d= (rposition_elec(:,4)-rposition_elec(:,3))/2;
0435
0436
0437 depth= depth+depthLevel;
0438 P= depth+1i*location;
0439 [Pu,iu,ju]= unique(round(P/depthPrecision)*depthPrecision);
0440
0441 if length(Pu) < length(P)
0442 lu= location(iu);
0443 zu= depth(iu);
0444 du= iu*0;
0445 for i= 1:length(Pu)
0446 du(i)= mean(data(ju==i));
0447 end
0448 if length(fmdl.misc.sizepoint)>length(Pu)
0449 su= iu*0;
0450
0451
0452 for i= 1:length(Pu)
0453 su(i)= min(fmdl.misc.sizepoint(ju==i));
0454
0455
0456 end
0457 fmdl.misc.sizepoint= su;
0458
0459 end
0460 else
0461 lu= location;
0462 zu= depth;
0463 du= data;
0464 end
0465
0466 [Puu]= unique(lu+1i*zu);
0467 if length(Puu) ~= length(lu)
0468 keyboard
0469 end
0470 scatter(lu,zu,fmdl.misc.sizepoint,(du),'filled','MarkerEdgeColor','k');
0471
0472 xlabel('Distance (m)','fontsize',fs,'fontname','Times','interpreter','latex');
0473 ylabel('Pseudo-depth (m)','fontsize',fs,'fontname','Times','interpreter','latex')
0474
0475 set(gca,'fontsize',fs,'fontname','Times')
0476 end
0477
0478
0479
0480 function [zps,xps]= plotPseudoSectionProfileDown(fmdl,data)
0481 fs= fmdl.show_pseudosection.fs;
0482
0483 [elec_posn,elecNumber]= electrodesPosition(fmdl);
0484 xposition_elec= reshape(elec_posn(elecNumber,1),[],4);
0485 xps= mean(xposition_elec,2);
0486
0487 AB= abs(xposition_elec(:,2)-xposition_elec(:,1));
0488 MN= abs(xposition_elec(:,4)-xposition_elec(:,3));
0489 AN= abs(xposition_elec(:,4)-xposition_elec(:,1));
0490 BM= abs(xposition_elec(:,3)-xposition_elec(:,2));
0491 AM= abs(xposition_elec(:,3)-xposition_elec(:,1));
0492 BN= abs(xposition_elec(:,4)-xposition_elec(:,2));
0493
0494 a= MN;
0495 zps= -a/3;
0496
0497
0498
0499
0500
0501
0502
0503
0504
0505
0506
0507
0508 P= zps+1i*xps;
0509 [Pu,iu,ju]= unique(P);
0510
0511 if length(Pu) < length(P)
0512 xu= xps(iu);
0513 zu= zps(iu);
0514 du= iu*0;
0515 for i= 1:length(Pu)
0516 du(i)= mean(data(ju==i));
0517 end
0518 else
0519 xu= xps;
0520 zu= zps;
0521 du= data;
0522 end
0523
0524 scatter(xu,zu,fmdl.misc.sizepoint,(du),'filled','MarkerEdgeColor','k');
0525 xlabel('Distance (m)','fontsize',fs,'fontname','Times','interpreter','latex');
0526 ylabel('Pseudo-depth (m)','fontsize',fs,'fontname','Times','interpreter','latex');
0527 axis equal; axis tight;
0528 bx= get(gca,'xlim'); by= get(gca,'ylim');
0529 set(gca,'fontsize',fs,'fontname','Times','dataAspectRatio',[1 (by(2)-by(1))/(bx(2)-bx(1))*3 1])
0530
0531 end
0532
0533 function [xps,zps]= plotPseudoSectionProfileVert(fmdl,data)
0534 fs= fmdl.show_pseudosection.fs;
0535
0536 [elec_posn,elecNumber]= electrodesPosition(fmdl);
0537
0538 zposition_elec= reshape(elec_posn(elecNumber,3),[],4);
0539 zps= mean(zposition_elec,2);
0540
0541 AB= abs(zposition_elec(:,2)-zposition_elec(:,1));
0542 MN= abs(zposition_elec(:,4)-zposition_elec(:,3));
0543 AN= abs(zposition_elec(:,4)-zposition_elec(:,1));
0544 BM= abs(zposition_elec(:,3)-zposition_elec(:,2));
0545 AM= abs(zposition_elec(:,3)-zposition_elec(:,1));
0546 BN= abs(zposition_elec(:,4)-zposition_elec(:,2));
0547
0548 a= max([AB MN AN BM AM BN ],[],2);
0549 xps= a/3;
0550
0551 P= zps+1i*xps;
0552 [Pu,iu,ju]= unique(P);
0553
0554 xu= xps(iu);
0555 zu= zps(iu);
0556 du= iu*0;
0557 for i= 1:length(Pu)
0558 du(i)= mean(data(ju==i));
0559 end
0560
0561 scatter(xu,zu,fmdl.misc.sizepoint,(du),'filled','MarkerEdgeColor','k');
0562 xlabel('Pseudo distance (m)','fontsize',fs,'fontname','Times','interpreter','latex');
0563 if zps(1)<0
0564 ylabel('Depth (m)','fontsize',fs,'fontname','Times','interpreter','latex')
0565 else
0566 ylabel('Height (m)','fontsize',fs,'fontname','Times','interpreter','latex')
0567 end
0568 axis equal; axis tight;
0569 set(gca,'fontsize',fs,'fontname','Times')
0570 end
0571
0572 function [r_point,th_point]= plotPseudoSectionCircularIn(fmdl,data)
0573 fs= fmdl.show_pseudosection.fs;
0574 [elec_posn,elecNumber] = electrodesPosition(fmdl);
0575
0576 [A,~,th_point,r,xc,yc] = polarPosition(elecNumber,elec_posn);
0577 a= max(A,[],2);
0578 r_point= a/2;
0579 r_point= (r-r_point/pi);
0580
0581 [x_point,y_point]= pol2cart(th_point,r_point);
0582 x_point= x_point+xc; y_point= y_point+yc;
0583
0584 P= x_point+1i*y_point;
0585 [Pu,iu,ju]= unique(P);
0586
0587 xu= x_point(iu);
0588 zu= y_point(iu);
0589 du= iu*0;
0590 for i= 1:length(Pu)
0591 du(i)= mean(data(ju==i));
0592 end
0593
0594 scatter(xu,zu,fmdl.misc.sizepoint,(du),'filled','MarkerEdgeColor','k');
0595 xlabel('X (m)','fontsize',fs,'fontname','Times','interpreter','latex');
0596 ylabel('Y (m)','fontsize',fs,'fontname','Times','interpreter','latex')
0597 axis equal; axis tight;
0598 set(gca,'fontsize',fs,'fontname','Times')
0599 end
0600
0601 function [dMN,th_point]= plotPseudoSectionCircularInVolcano(fmdl,data)
0602 fs= fmdl.show_pseudosection.fs;
0603 [elec_posn,elecNumber] = electrodesPosition(fmdl);
0604 [xposition_elec,yposition_elec,zposition_elec] = electrodesPositionABMN(elecNumber,elec_posn);
0605 dMN= sqrt((xposition_elec(:,4)-xposition_elec(:,3)).^2 + ...
0606 (yposition_elec(:,4)-yposition_elec(:,3)).^2 + (zposition_elec(:,4)-zposition_elec(:,3)).^2);
0607
0608 [A,r_bary,th_point,r,xc,yc,TH,R] = polarPosition(elecNumber,elec_posn);
0609 [th_bary,r_bary,th_bary34,r_bary34,r,xc,yc,barycenter_x34,barycenter_y34] = bary34(elecNumber,elec_posn);
0610
0611
0612 r= 400;
0613 r_point= (r-r_bary)/(max(r-r_bary)-min(r-r_bary))*(r-5);
0614 r_point= r_point-min(r_point)+5;
0615 [x_point,y_point]= pol2cart(th_bary34,r_point);
0616 x_point= x_point+xc; y_point= y_point+yc;
0617 P= x_point+1i*y_point;
0618 [Pu,iu,ju]= unique(P);
0619 xu= x_point(iu); zu= y_point(iu); du= iu*0;
0620 for i= 1:length(Pu); du(i)= mean(data(ju==i)); end
0621 figure; scatter(xu,zu,fmdl.misc.sizepoint,(du),'filled','MarkerEdgeColor','k');
0622 text(xu(iu==1),zu(iu==1),'1')
0623 hold on; plot(xposition_elec,yposition_elec,'kp')
0624 xlabel('X (m)','fontsize',fs,'fontname','Times','interpreter','latex');
0625 ylabel('Y (m)','fontsize',fs,'fontname','Times','interpreter','latex')
0626 axis equal; axis tight;
0627 set(gca,'fontsize',fs,'fontname','Times')
0628 end
0629
0630 function [r_point,th_point]= plotPseudoSectionCircularOut(fmdl,data)
0631 fs= fmdl.show_pseudosection.fs;
0632 [elec_posn,elecNumber] = electrodesPosition(fmdl);
0633 [A,r_bary,th_point,r,xc,yc] = polarPosition(elecNumber,elec_posn);
0634 a= max(A,[],2);
0635 r_point= a/2;
0636 r_point= (r+r_point);
0637 [x_point,y_point]= pol2cart(th_point,r_point+r);
0638 x_point= x_point+xc; y_point= y_point+yc;
0639
0640 P= x_point+1i*y_point;
0641 [Pu,iu,ju]= unique(P);
0642
0643 xu= x_point(iu);
0644 zu= y_point(iu);
0645 du= iu*0;
0646 for i= 1:length(Pu)
0647 du(i)= mean(data(ju==i));
0648 end
0649
0650 scatter(xu,zu,fmdl.misc.sizepoint,(du),'filled','MarkerEdgeColor','k'); colorbar
0651 xlabel('X (m)','fontsize',fs,'fontname','Times','interpreter','latex');
0652 ylabel('Y (m)','fontsize',fs,'fontname','Times','interpreter','latex')
0653 axis equal; axis tight;
0654 set(gca,'fontsize',fs,'fontname','Times')
0655 end
0656
0657
0658 function [elec_posn,elecNumber] = electrodesPosition(fmdl)
0659 stim= fmdl.stimulation;
0660 stimulationMatrix= [];
0661 for i = 1:length(stim);
0662 nmp= size(stim(i).meas_pattern, 1);
0663 [idxIN,idxJN]= find(stim(i).meas_pattern<0);
0664 [idxIP,idxJP]= find(stim(i).meas_pattern>0);
0665 stimulationMatrix= [stimulationMatrix; [ 0*(1:nmp)'+find(stim(i).stim_pattern<0) ...
0666 0*(1:nmp)'+find(stim(i).stim_pattern>0)] idxJN(idxIN) idxJP(idxIP)];
0667 end
0668
0669
0670
0671
0672
0673
0674
0675
0676
0677
0678
0679
0680
0681
0682
0683
0684
0685 elecNumber= stimulationMatrix;
0686
0687 elec_posn= zeros(length(fmdl.electrode),size(fmdl.nodes,2));
0688 for i=1:length(fmdl.electrode)
0689 elec_posn(i,:)= mean(fmdl.nodes(fmdl.electrode(1,i).nodes,:),1);
0690 end
0691
0692 end
0693
0694 function [xposition_elec,yposition_elec,zposition_elec] = electrodesPositionABMN(elecNumber,elec_posn)
0695 xposition_elec= reshape(elec_posn(elecNumber,1),[],4);
0696 yposition_elec= reshape(elec_posn(elecNumber,2),[],4);
0697 zposition_elec= reshape(elec_posn(elecNumber,3),[],4);
0698 end
0699
0700 function [A,r_bary,th_bary,r,xc,yc,TH,R] = polarPosition(elecNumber,elec_posn)
0701 TH= elecNumber*0;
0702 R= elecNumber*0;
0703
0704 xposition_elec= reshape(elec_posn(elecNumber,1),[],4);
0705 yposition_elec= reshape(elec_posn(elecNumber,2),[],4);
0706 rx= (max(xposition_elec(:))-min(xposition_elec(:)))/2;
0707 ry= (max(yposition_elec(:))-min(yposition_elec(:)))/2;
0708 rmean= (rx+ry)/2;
0709
0710 xc= mean(elec_posn(:,1)); yc= mean(elec_posn(:,2)); r = rmean;
0711
0712 barycenter_x= mean(xposition_elec,2);
0713 barycenter_y= mean(yposition_elec,2);
0714
0715 [th_bary,r_bary]= cart2pol(barycenter_x-xc,barycenter_y-yc);
0716
0717 for i= 1:4
0718 [TH(:,i),R(:,i)]= cart2pol(elec_posn(elecNumber(:,i),1)-xc,elec_posn(elecNumber(:,i),2)-yc);
0719 end
0720
0721 Arc_length_AB= (TH(:,2)-TH(:,1));
0722 idx_discontinuity= find(TH(:,2)>TH(:,1) & TH(:,2)<0);
0723 Arc_length_AB(idx_discontinuity)= (-TH(idx_discontinuity,1)+TH(idx_discontinuity,2));
0724 idx_discontinuity= find(TH(:,1)>TH(:,2));
0725 Arc_length_AB(idx_discontinuity)= (TH(idx_discontinuity,1)-TH(idx_discontinuity,2));
0726 idx_discontinuity= find(TH(:,1)>TH(:,2) & TH(:,1)<0);
0727 Arc_length_AB(idx_discontinuity)= (-TH(idx_discontinuity,2)+TH(idx_discontinuity,1));
0728 Arc_length_AB(Arc_length_AB>=pi+0.005)= 2*pi-Arc_length_AB(Arc_length_AB>=pi+0.005);
0729 Arc_length_AB= Arc_length_AB.*(R(:,1)+R(:,2))/2;
0730
0731 Arc_length_MN= (TH(:,4)-TH(:,3));
0732 idx_discontinuity= find((TH(:,4)>TH(:,3) & TH(:,4)<0));
0733 Arc_length_MN(idx_discontinuity)= -TH(idx_discontinuity,3)+TH(idx_discontinuity,4);
0734 idx_discontinuity= find((TH(:,4)<TH(:,3) & TH(:,4)<0));
0735 Arc_length_MN(idx_discontinuity)= +TH(idx_discontinuity,3)-TH(idx_discontinuity,4);
0736 idx_discontinuity= find((TH(:,3)>TH(:,4)) & TH(:,3)>0);
0737 Arc_length_MN(idx_discontinuity)= TH(idx_discontinuity,3)-TH(idx_discontinuity,4);
0738
0739
0740 Arc_length_MN(Arc_length_MN>=pi+0.005)= 2*pi-Arc_length_MN(Arc_length_MN>=pi+0.005);
0741 Arc_length_MN= Arc_length_MN.*(R(:,3)+R(:,4))/2;
0742
0743 Arc_length_AM= (TH(:,3)-TH(:,1));
0744 idx_discontinuity= find((TH(:,3)>TH(:,1) & TH(:,3)<0));
0745 Arc_length_AM(idx_discontinuity)= -TH(idx_discontinuity,1)+TH(idx_discontinuity,3);
0746 idx_discontinuity= find((TH(:,1)> TH(:,3)));
0747 Arc_length_AM(idx_discontinuity)= TH(idx_discontinuity,1)-TH(idx_discontinuity,3);
0748 idx_discontinuity= find((TH(:,1)> TH(:,3) & TH(:,1) <0));
0749 Arc_length_AM(idx_discontinuity)= -TH(idx_discontinuity,1)+TH(idx_discontinuity,3);
0750 Arc_length_AM(Arc_length_AM>=pi+0.005)= 2*pi-Arc_length_AM(Arc_length_AM>=pi+0.005);
0751 Arc_length_AM= Arc_length_AM.*(R(:,1)+R(:,3))/2;
0752
0753 Arc_length_AN= (TH(:,4)-TH(:,1));
0754 idx_discontinuity= find((TH(:,4)>TH(:,1)<0 & TH(:,4)<0));
0755 Arc_length_AN(idx_discontinuity)= -TH(idx_discontinuity,1)+TH(idx_discontinuity,4);
0756 idx_discontinuity= find((TH(:,1)>TH(:,4)));
0757 Arc_length_AN(idx_discontinuity)= TH(idx_discontinuity,1)-TH(idx_discontinuity,4);
0758 idx_discontinuity= find((TH(:,1)>TH(:,4) & TH(:,1)<0));
0759 Arc_length_AN(idx_discontinuity)= TH(idx_discontinuity,4)-TH(idx_discontinuity,1);
0760 Arc_length_AN(Arc_length_AN>=pi+0.005)= 2*pi-Arc_length_AN(Arc_length_AN>=pi+0.005);
0761 Arc_length_AN= Arc_length_AN.*(R(:,1)+R(:,4))/2;
0762
0763 Arc_length_NB= (TH(:,4)-TH(:,2));
0764 idx_discontinuity= find((TH(:,4)>TH(:,2)));
0765 Arc_length_NB(idx_discontinuity)= -TH(idx_discontinuity,2)+TH(idx_discontinuity,4);
0766 idx_discontinuity= find((TH(:,2)> TH(:,4)));
0767 Arc_length_NB(idx_discontinuity)= -TH(idx_discontinuity,4)+TH(idx_discontinuity,2);
0768 idx_discontinuity= find((TH(:,2)> TH(:,4) & TH(:,4)<0));
0769 Arc_length_NB(idx_discontinuity)= -TH(idx_discontinuity,4)+TH(idx_discontinuity,2);
0770 Arc_length_NB(Arc_length_NB>=pi+0.005)= 2*pi-Arc_length_NB(Arc_length_NB>=pi+0.005);
0771 Arc_length_NB= Arc_length_NB.*(R(:,2)+R(:,4))/2;
0772
0773 Arc_length_BM= (TH(:,3)-TH(:,2));
0774 idx_discontinuity= find((TH(:,3)>TH(:,2) & TH(:,2)<0));
0775 Arc_length_BM(idx_discontinuity)= -TH(idx_discontinuity,2)+TH(idx_discontinuity,3);
0776 idx_discontinuity= find((TH(:,2)> TH(:,3)));
0777 Arc_length_BM(idx_discontinuity)= -TH(idx_discontinuity,3)+TH(idx_discontinuity,2);
0778 idx_discontinuity= find((TH(:,2)> TH(:,3)& TH(:,3)<0));
0779 Arc_length_BM(idx_discontinuity)= TH(idx_discontinuity,3)-TH(idx_discontinuity,2);
0780 Arc_length_BM(Arc_length_BM>=pi+0.005)= 2*pi-Arc_length_BM(Arc_length_BM>=pi+0.005);
0781 Arc_length_BM= Arc_length_BM.*(R(:,2)+R(:,3))/2;
0782
0783 A= [Arc_length_AB Arc_length_MN Arc_length_AM Arc_length_AN Arc_length_NB Arc_length_BM];
0784 end
0785
0786
0787
0788 function [th_bary,r_bary,th_bary34,r_bary34,r,xc,yc,barycenter_x34,barycenter_y34] = bary34(elecNumber,elec_posn)
0789 xposition_elec= reshape(elec_posn(elecNumber,1),[],4);
0790 yposition_elec= reshape(elec_posn(elecNumber,2),[],4);
0791 rx= (max(xposition_elec(:))-min(xposition_elec(:)))/2;
0792 ry= (max(yposition_elec(:))-min(yposition_elec(:)))/2;
0793 rmean= (rx+ry)/2;
0794
0795 xc= mean(elec_posn(:,1)); yc= mean(elec_posn(:,2)); r = rmean;
0796
0797 barycenter_x34= mean(xposition_elec(:,3:4),2);
0798 barycenter_y34= mean(yposition_elec(:,3:4),2);
0799 [th_bary34,r_bary34]= cart2pol(barycenter_x34-xc,barycenter_y34-yc);
0800
0801 barycenter_x= mean(xposition_elec,2);
0802 barycenter_y= mean(yposition_elec,2);
0803 [th_bary,r_bary]= cart2pol(barycenter_x-xc,barycenter_y-yc);
0804
0805
0806 end
0807
0808 function do_unit_test
0809 shape_str = ['solid top = plane(0,0,0;0,1,0);\n' ...
0810 'solid mainobj= top and orthobrick(-100,-200,-100;410,10,100) -maxh=20.0;\n'];
0811 e0 = linspace(0,310,64)';
0812 elec_pos = [e0,0*e0,0*e0,1+0*e0,0*e0,0*e0];
0813 elec_shape= [0.1,0.1,1];
0814 elec_obj = 'top';
0815 fmdl = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0816
0817 spacing= [1 1 1 2 3 3 4 4 5 6 6 7 8 8 9 10 10 11 12 12 13 14 14 15 16 17];
0818 multiples= [1 2 3 2 1 5/3 1 2 1 1 7/6 1 1 10/8 1 1 12/10 1 1 13/12 1 1 15/14 1 1 1];
0819 fmdl.stimulation= stim_pattern_geophys( 64, 'DipoleDipole', {'spacings', spacing,'multiples',multiples} );
0820
0821 img1= mk_image(fmdl,1);
0822 vh1= fwd_solve(img1);
0823 normalisation= 1./vh1.meas;
0824 I= speye(length(normalisation));
0825 I(1:size(I,1)+1:size(I,1)*size(I,1))= normalisation;
0826
0827 img = mk_image(fmdl,0+ mk_c2f_circ_mapping(fmdl,[100;-50;0;50])*100);
0828 img.elem_data(img.elem_data==0)= 0.1;
0829 dd = fwd_solve(img);
0830 subplot(221);
0831 show_pseudosection( fmdl, I*dd.meas )
0832
0833 subplot(222);
0834 fmdl.show_pseudosection.orientation = 'horizontaldownward';
0835 show_pseudosection( fmdl, I*dd.meas )
0836
0837 subplot(223);
0838 fmdl.show_pseudosection.orientation = 'vertical';
0839 show_pseudosection( fmdl, I*dd.meas )
0840
0841 end