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 [A,th_point,r,xc,yc] = polarPosition(elecNumber,elec_posn);
0576 a= max(A,[],2);
0577 r_point= a/2;
0578 r_point= (r-r_point/pi);
0579
0580 [x_point,y_point]= pol2cart(th_point,r_point);
0581 x_point= x_point+xc; y_point= y_point+yc;
0582
0583 P= x_point+1i*y_point;
0584 [Pu,iu,ju]= unique(P);
0585
0586 xu= x_point(iu);
0587 zu= y_point(iu);
0588 du= iu*0;
0589 for i= 1:length(Pu)
0590 du(i)= mean(data(ju==i));
0591 end
0592
0593 scatter(xu,zu,fmdl.misc.sizepoint,(du),'filled','MarkerEdgeColor','k');
0594 xlabel('X (m)','fontsize',fs,'fontname','Times','interpreter','latex');
0595 ylabel('Y (m)','fontsize',fs,'fontname','Times','interpreter','latex')
0596 axis equal; axis tight;
0597 set(gca,'fontsize',fs,'fontname','Times')
0598 end
0599
0600 function [dMN,th_point]= plotPseudoSectionCircularInVolcano(fmdl,data)
0601 fs= fmdl.show_pseudosection.fs;
0602 [elec_posn,elecNumber] = electrodesPosition(fmdl);
0603 [xposition_elec,yposition_elec,zposition_elec] = electrodesPositionABMN(elecNumber,elec_posn);
0604 dMN= sqrt((xposition_elec(:,4)-xposition_elec(:,3)).^2 + ...
0605 (yposition_elec(:,4)-yposition_elec(:,3)).^2 + (zposition_elec(:,4)-zposition_elec(:,3)).^2);
0606
0607 [A,r_bary,th_point,r,xc,yc,TH,R] = polarPosition(elecNumber,elec_posn);
0608 [th_bary,r_bary,th_bary34,r_bary34,r,xc,yc,barycenter_x34,barycenter_y34] = bary34(elecNumber,elec_posn);
0609
0610
0611 r= 400;
0612 r_point= (r-r_bary)/(max(r-r_bary)-min(r-r_bary))*(r-5);
0613 r_point= r_point-min(r_point)+5;
0614 [x_point,y_point]= pol2cart(th_bary34,r_point);
0615 x_point= x_point+xc; y_point= y_point+yc;
0616 P= x_point+1i*y_point;
0617 [Pu,iu,ju]= unique(P);
0618 xu= x_point(iu); zu= y_point(iu); du= iu*0;
0619 for i= 1:length(Pu); du(i)= mean(data(ju==i)); end
0620 figure; scatter(xu,zu,fmdl.misc.sizepoint,(du),'filled','MarkerEdgeColor','k');
0621 text(xu(iu==1),zu(iu==1),'1')
0622 hold on; plot(xposition_elec,yposition_elec,'kp')
0623 xlabel('X (m)','fontsize',fs,'fontname','Times','interpreter','latex');
0624 ylabel('Y (m)','fontsize',fs,'fontname','Times','interpreter','latex')
0625 axis equal; axis tight;
0626 set(gca,'fontsize',fs,'fontname','Times')
0627 end
0628
0629 function [r_point,th_point]= plotPseudoSectionCircularOut(fmdl,data)
0630 fs= fmdl.show_pseudosection.fs;
0631 [elec_posn,elecNumber] = electrodesPosition(fmdl);
0632 [A,r_bary,th_point,r,xc,yc] = polarPosition(elecNumber,elec_posn);
0633 a= max(A,[],2);
0634 r_point= a/2;
0635 r_point= (r+r_point);
0636 [x_point,y_point]= pol2cart(th_point,r_point+r);
0637 x_point= x_point+xc; y_point= y_point+yc;
0638
0639 P= x_point+1i*y_point;
0640 [Pu,iu,ju]= unique(P);
0641
0642 xu= x_point(iu);
0643 zu= y_point(iu);
0644 du= iu*0;
0645 for i= 1:length(Pu)
0646 du(i)= mean(data(ju==i));
0647 end
0648
0649 scatter(xu,zu,fmdl.misc.sizepoint,(du),'filled','MarkerEdgeColor','k'); colorbar
0650 xlabel('X (m)','fontsize',fs,'fontname','Times','interpreter','latex');
0651 ylabel('Y (m)','fontsize',fs,'fontname','Times','interpreter','latex')
0652 axis equal; axis tight;
0653 set(gca,'fontsize',fs,'fontname','Times')
0654 end
0655
0656
0657 function [elec_posn,elecNumber] = electrodesPosition(fmdl)
0658 stim= fmdl.stimulation;
0659 stimulationMatrix= [];
0660 for i = 1:length(stim);
0661 nmp= size(stim(i).meas_pattern, 1);
0662 [idxIN,idxJN]= find(stim(i).meas_pattern<0);
0663 [idxIP,idxJP]= find(stim(i).meas_pattern>0);
0664 stimulationMatrix= [stimulationMatrix; [ 0*(1:nmp)'+find(stim(i).stim_pattern<0) ...
0665 0*(1:nmp)'+find(stim(i).stim_pattern>0)] idxJN(idxIN) idxJP(idxIP)];
0666 end
0667
0668
0669
0670
0671
0672
0673
0674
0675
0676
0677
0678
0679
0680
0681
0682
0683
0684 elecNumber= stimulationMatrix;
0685
0686 elec_posn= zeros(length(fmdl.electrode),size(fmdl.nodes,2));
0687 for i=1:length(fmdl.electrode)
0688 elec_posn(i,:)= mean(fmdl.nodes(fmdl.electrode(1,i).nodes,:),1);
0689 end
0690
0691 end
0692
0693 function [xposition_elec,yposition_elec,zposition_elec] = electrodesPositionABMN(elecNumber,elec_posn)
0694 xposition_elec= reshape(elec_posn(elecNumber,1),[],4);
0695 yposition_elec= reshape(elec_posn(elecNumber,2),[],4);
0696 zposition_elec= reshape(elec_posn(elecNumber,3),[],4);
0697 end
0698
0699 function [A,r_bary,th_bary,r,xc,yc,TH,R] = polarPosition(elecNumber,elec_posn)
0700 TH= elecNumber*0;
0701 R= elecNumber*0;
0702
0703 xposition_elec= reshape(elec_posn(elecNumber,1),[],4);
0704 yposition_elec= reshape(elec_posn(elecNumber,2),[],4);
0705 rx= (max(xposition_elec(:))-min(xposition_elec(:)))/2;
0706 ry= (max(yposition_elec(:))-min(yposition_elec(:)))/2;
0707 rmean= (rx+ry)/2;
0708
0709 xc= mean(elec_posn(:,1)); yc= mean(elec_posn(:,2)); r = rmean;
0710
0711 barycenter_x= mean(xposition_elec,2);
0712 barycenter_y= mean(yposition_elec,2);
0713
0714 [th_bary,r_bary]= cart2pol(barycenter_x-xc,barycenter_y-yc);
0715
0716 for i= 1:4
0717 [TH(:,i),R(:,i)]= cart2pol(elec_posn(elecNumber(:,i),1)-xc,elec_posn(elecNumber(:,i),2)-yc);
0718 end
0719
0720 Arc_length_AB= (TH(:,2)-TH(:,1));
0721 idx_discontinuity= find(TH(:,2)>TH(:,1) & TH(:,2)<0);
0722 Arc_length_AB(idx_discontinuity)= (-TH(idx_discontinuity,1)+TH(idx_discontinuity,2));
0723 idx_discontinuity= find(TH(:,1)>TH(:,2));
0724 Arc_length_AB(idx_discontinuity)= (TH(idx_discontinuity,1)-TH(idx_discontinuity,2));
0725 idx_discontinuity= find(TH(:,1)>TH(:,2) & TH(:,1)<0);
0726 Arc_length_AB(idx_discontinuity)= (-TH(idx_discontinuity,2)+TH(idx_discontinuity,1));
0727 Arc_length_AB(Arc_length_AB>=pi+0.005)= 2*pi-Arc_length_AB(Arc_length_AB>=pi+0.005);
0728 Arc_length_AB= Arc_length_AB.*(R(:,1)+R(:,2))/2;
0729
0730 Arc_length_MN= (TH(:,4)-TH(:,3));
0731 idx_discontinuity= find((TH(:,4)>TH(:,3) & TH(:,4)<0));
0732 Arc_length_MN(idx_discontinuity)= -TH(idx_discontinuity,3)+TH(idx_discontinuity,4);
0733 idx_discontinuity= find((TH(:,4)<TH(:,3) & TH(:,4)<0));
0734 Arc_length_MN(idx_discontinuity)= +TH(idx_discontinuity,3)-TH(idx_discontinuity,4);
0735 idx_discontinuity= find((TH(:,3)>TH(:,4)) & TH(:,3)>0);
0736 Arc_length_MN(idx_discontinuity)= TH(idx_discontinuity,3)-TH(idx_discontinuity,4);
0737
0738
0739 Arc_length_MN(Arc_length_MN>=pi+0.005)= 2*pi-Arc_length_MN(Arc_length_MN>=pi+0.005);
0740 Arc_length_MN= Arc_length_MN.*(R(:,3)+R(:,4))/2;
0741
0742 Arc_length_AM= (TH(:,3)-TH(:,1));
0743 idx_discontinuity= find((TH(:,3)>TH(:,1) & TH(:,3)<0));
0744 Arc_length_AM(idx_discontinuity)= -TH(idx_discontinuity,1)+TH(idx_discontinuity,3);
0745 idx_discontinuity= find((TH(:,1)> TH(:,3)));
0746 Arc_length_AM(idx_discontinuity)= TH(idx_discontinuity,1)-TH(idx_discontinuity,3);
0747 idx_discontinuity= find((TH(:,1)> TH(:,3) & TH(:,1) <0));
0748 Arc_length_AM(idx_discontinuity)= -TH(idx_discontinuity,1)+TH(idx_discontinuity,3);
0749 Arc_length_AM(Arc_length_AM>=pi+0.005)= 2*pi-Arc_length_AM(Arc_length_AM>=pi+0.005);
0750 Arc_length_AM= Arc_length_AM.*(R(:,1)+R(:,3))/2;
0751
0752 Arc_length_AN= (TH(:,4)-TH(:,1));
0753 idx_discontinuity= find((TH(:,4)>TH(:,1)<0 & TH(:,4)<0));
0754 Arc_length_AN(idx_discontinuity)= -TH(idx_discontinuity,1)+TH(idx_discontinuity,4);
0755 idx_discontinuity= find((TH(:,1)>TH(:,4)));
0756 Arc_length_AN(idx_discontinuity)= TH(idx_discontinuity,1)-TH(idx_discontinuity,4);
0757 idx_discontinuity= find((TH(:,1)>TH(:,4) & TH(:,1)<0));
0758 Arc_length_AN(idx_discontinuity)= TH(idx_discontinuity,4)-TH(idx_discontinuity,1);
0759 Arc_length_AN(Arc_length_AN>=pi+0.005)= 2*pi-Arc_length_AN(Arc_length_AN>=pi+0.005);
0760 Arc_length_AN= Arc_length_AN.*(R(:,1)+R(:,4))/2;
0761
0762 Arc_length_NB= (TH(:,4)-TH(:,2));
0763 idx_discontinuity= find((TH(:,4)>TH(:,2)));
0764 Arc_length_NB(idx_discontinuity)= -TH(idx_discontinuity,2)+TH(idx_discontinuity,4);
0765 idx_discontinuity= find((TH(:,2)> TH(:,4)));
0766 Arc_length_NB(idx_discontinuity)= -TH(idx_discontinuity,4)+TH(idx_discontinuity,2);
0767 idx_discontinuity= find((TH(:,2)> TH(:,4) & TH(:,4)<0));
0768 Arc_length_NB(idx_discontinuity)= -TH(idx_discontinuity,4)+TH(idx_discontinuity,2);
0769 Arc_length_NB(Arc_length_NB>=pi+0.005)= 2*pi-Arc_length_NB(Arc_length_NB>=pi+0.005);
0770 Arc_length_NB= Arc_length_NB.*(R(:,2)+R(:,4))/2;
0771
0772 Arc_length_BM= (TH(:,3)-TH(:,2));
0773 idx_discontinuity= find((TH(:,3)>TH(:,2) & TH(:,2)<0));
0774 Arc_length_BM(idx_discontinuity)= -TH(idx_discontinuity,2)+TH(idx_discontinuity,3);
0775 idx_discontinuity= find((TH(:,2)> TH(:,3)));
0776 Arc_length_BM(idx_discontinuity)= -TH(idx_discontinuity,3)+TH(idx_discontinuity,2);
0777 idx_discontinuity= find((TH(:,2)> TH(:,3)& TH(:,3)<0));
0778 Arc_length_BM(idx_discontinuity)= TH(idx_discontinuity,3)-TH(idx_discontinuity,2);
0779 Arc_length_BM(Arc_length_BM>=pi+0.005)= 2*pi-Arc_length_BM(Arc_length_BM>=pi+0.005);
0780 Arc_length_BM= Arc_length_BM.*(R(:,2)+R(:,3))/2;
0781
0782 A= [Arc_length_AB Arc_length_MN Arc_length_AM Arc_length_AN Arc_length_NB Arc_length_BM];
0783 end
0784
0785
0786
0787 function [th_bary,r_bary,th_bary34,r_bary34,r,xc,yc,barycenter_x34,barycenter_y34] = bary34(elecNumber,elec_posn)
0788 xposition_elec= reshape(elec_posn(elecNumber,1),[],4);
0789 yposition_elec= reshape(elec_posn(elecNumber,2),[],4);
0790 rx= (max(xposition_elec(:))-min(xposition_elec(:)))/2;
0791 ry= (max(yposition_elec(:))-min(yposition_elec(:)))/2;
0792 rmean= (rx+ry)/2;
0793
0794 xc= mean(elec_posn(:,1)); yc= mean(elec_posn(:,2)); r = rmean;
0795
0796 barycenter_x34= mean(xposition_elec(:,3:4),2);
0797 barycenter_y34= mean(yposition_elec(:,3:4),2);
0798 [th_bary34,r_bary34]= cart2pol(barycenter_x34-xc,barycenter_y34-yc);
0799
0800 barycenter_x= mean(xposition_elec,2);
0801 barycenter_y= mean(yposition_elec,2);
0802 [th_bary,r_bary]= cart2pol(barycenter_x-xc,barycenter_y-yc);
0803
0804
0805 end
0806
0807 function do_unit_test
0808 shape_str = ['solid top = plane(0,0,0;0,1,0);\n' ...
0809 'solid mainobj= top and orthobrick(-100,-200,-100;410,10,100) -maxh=20.0;\n'];
0810 e0 = linspace(0,310,64)';
0811 elec_pos = [e0,0*e0,0*e0,1+0*e0,0*e0,0*e0];
0812 elec_shape= [0.1,0.1,1];
0813 elec_obj = 'top';
0814 fmdl = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0815
0816 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];
0817 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];
0818 fmdl.stimulation= stim_pattern_geophys( 64, 'DipoleDipole', {'spacings', spacing,'multiples',multiples} );
0819
0820 img1= mk_image(fmdl,1);
0821 vh1= fwd_solve(img1);
0822 normalisation= 1./vh1.meas;
0823 I= speye(length(normalisation));
0824 I(1:size(I,1)+1:size(I,1)*size(I,1))= normalisation;
0825
0826 img = mk_image(fmdl,0+ mk_c2f_circ_mapping(fmdl,[100;-50;0;50])*100);
0827 img.elem_data(img.elem_data==0)= 0.1;
0828 dd = fwd_solve(img);
0829 subplot(221);
0830 show_pseudosection( fmdl, I*dd.meas )
0831
0832 subplot(222);
0833 fmdl.show_pseudosection.orientation = 'horizontaldownward';
0834 show_pseudosection( fmdl, I*dd.meas )
0835
0836 subplot(223);
0837 fmdl.show_pseudosection.orientation = 'vertical';
0838 show_pseudosection( fmdl, I*dd.meas )
0839
0840 end