show_pseudosection

PURPOSE ^

SHOW_PSEUDOSECTION: show a pseudo-section image of data

SYNOPSIS ^

function fwd_model= show_pseudosection( fwd_model, data, orientation)

DESCRIPTION ^

SHOW_PSEUDOSECTION: show a pseudo-section image of data

 OPTIONS:
   fwd_model - fwd model of the domain
   data      - data to display

 OPTIONS:
 fwd_model.show_pseudosection.orientation - of the pseudo-section
      orientation = 'horizontaldownward': section in +z direction from 
         the electrode plane (default)
      orientation = 'circularoutside': gallery outside electrode plane
      orientation = 'vertical'
      orientation = 'circularinside'

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function fwd_model= show_pseudosection( fwd_model, data, orientation)
0002 %SHOW_PSEUDOSECTION: show a pseudo-section image of data
0003 %
0004 % OPTIONS:
0005 %   fwd_model - fwd model of the domain
0006 %   data      - data to display
0007 %
0008 % OPTIONS:
0009 % fwd_model.show_pseudosection.orientation - of the pseudo-section
0010 %      orientation = 'horizontaldownward': section in +z direction from
0011 %         the electrode plane (default)
0012 %      orientation = 'circularoutside': gallery outside electrode plane
0013 %      orientation = 'vertical'
0014 %      orientation = 'circularinside'
0015 
0016 % (C) 2005-2008 Nolwenn Lesparre. License: GPL version 2 or version 3
0017 % $Id: show_pseudosection.m 4975 2015-05-11 01:53:02Z aadler $
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 % otherwise orientation is specified
0032 try
0033    orientation= fwd_model.show_pseudosection.orientation;
0034 catch
0035    orientation = 'horizontaldownward';
0036 end
0037 end
0038 
0039 if iscell(orientation) 
0040     if strcmp(orientation{2},'yz')
0041     fwd_model.nodes= fwd_model.nodes(:,[2 3 1]);
0042     rotation= 'yz';
0043     elseif strcmp(orientation{2},'xz')
0044      fwd_model.nodes= fwd_model.nodes(:,[1 3 2]);
0045      rotation= 'xz';
0046     end
0047     orientation= orientation{1};
0048 end
0049 
0050 switch(upper(orientation))
0051     case 'HORIZONTALDOWNWARD';  [depth,location]= plotPseudoSectionProfileDown(fwd_model,data);
0052     case 'VERTICAL';            [depth,location]= plotPseudoSectionProfileVert(fwd_model,data);
0053     case 'SPECIFICHORIZONTALDOWNWARD';  [depth,location]= plotPseudoSectionProfileSpecificDown(fwd_model,data); 
0054     case 'SPECIFICHORIZONTALUPWARD';  [depth,location,L,d]= plotPseudoSectionProfileSpecificUp(fwd_model,data);
0055     case 'CIRCULAROUTSIDE';     [depth,location]= plotPseudoSectionCircularOut(fwd_model,data);
0056     case 'CIRCULARINSIDE';      [depth,location]= plotPseudoSectionCircularIn(fwd_model,data);
0057     case 'CIRCULARVOLCANO';     [depth,location]= plotPseudoSectionCircularInVolcano(fwd_model,data);
0058   otherwise;
0059     error('No orientation of type "%s" available', upper(orientation));
0060 end
0061 
0062 if exist('rotation','var') 
0063     if strcmp(rotation,'yz')
0064     fwd_model.nodes= fwd_model.nodes(:,[3 1 2]);  
0065     elseif strcmp(rotation,'xz')
0066      fwd_model.nodes= fwd_model.nodes(:,[1 2 3]);
0067     end
0068 end
0069 
0070 fwd_model.show_pseudosection.depth= depth;
0071 fwd_model.show_pseudosection.location= location;
0072 % fwd_model.show_pseudosection.L= L;
0073 % fwd_model.show_pseudosection.d= d;
0074 end
0075 
0076 function [depth,location,electrodesLocation]= plotPseudoSectionProfileSpecificDown(fmdl,data)
0077    fs= 20;
0078    [elec_posn,elecNumber]= electrodesPosition(fmdl);
0079    xposition_elec= reshape(elec_posn(elecNumber,1),[],4); 
0080    yposition_elec= reshape(elec_posn(elecNumber,2),[],4); 
0081    zposition_elec= reshape(elec_posn(elecNumber,3),[],4); zposition_elec= zposition_elec-min(zposition_elec(:));
0082    if isfield(fmdl.show_pseudosection,'minx') && isfield(fmdl.show_pseudosection,'maxx') && ...
0083        isfield(fmdl.show_pseudosection,'miny') && isfield(fmdl.show_pseudosection,'maxy')
0084        minx= fmdl.show_pseudosection.minx;
0085        maxx= fmdl.show_pseudosection.maxx;
0086        miny= fmdl.show_pseudosection.miny;
0087        maxy= fmdl.show_pseudosection.maxy;
0088    else
0089        xposition_elec= xposition_elec-min(xposition_elec(:));
0090        yposition_elec= yposition_elec-min(yposition_elec(:));
0091        minx= min(xposition_elec(:));
0092        maxx= max(xposition_elec(:));
0093        miny= min(yposition_elec(:));
0094        maxy= max(yposition_elec(:));
0095    end
0096    
0097    if isfield(fmdl.show_pseudosection,'xToDeduce') && isfield(fmdl.show_pseudosection,'yToDeduce')
0098        if strcmp(fmdl.show_pseudosection.xToDeduce,'minx')
0099            xToDeduce= minx;
0100        elseif strcmp(fmdl.show_pseudosection.xToDeduce,'maxx')
0101            xToDeduce= maxx;
0102        end
0103        if strcmp(fmdl.show_pseudosection.yToDeduce,'miny')
0104            yToDeduce= miny;
0105        elseif strcmp(fmdl.show_pseudosection.yToDeduce,'maxy')
0106            yToDeduce= maxy;
0107        end
0108    else
0109        xToDeduce= minx;
0110        yToDeduce= miny;
0111    end
0112    rposition_elec= sqrt((xposition_elec-xToDeduce).^2 + ...
0113        (yposition_elec-yToDeduce).^2);
0114      
0115 %    rposition_elec= sqrt((xposition_elec-min(xposition_elec(:))).^2 + ...
0116 %        (yposition_elec-max(yposition_elec(:))).^2);
0117    
0118    if isfield(fmdl.show_pseudosection,'elecsUsed')
0119        elecsUsed= fmdl.show_pseudosection.elecsUsed;
0120    else
0121        elecsUsed= 3:4;
0122    end
0123    
0124    if isfield(fmdl.show_pseudosection,'dirAxis')
0125        dirAxis= fmdl.show_pseudosection.dirAxis;
0126    else
0127        dirAxis= 'X';
0128    end
0129    
0130    if isfield(fmdl.show_pseudosection,'depthRatio')
0131        depthRatio= fmdl.show_pseudosection.depthRatio;
0132    else
0133        depthRatio= 3;
0134    end
0135    
0136     if isfield(fmdl.show_pseudosection,'depthPrecision')
0137        depthPrecision= fmdl.show_pseudosection.depthPrecision;
0138    else
0139        depthPrecision= 1;
0140     end
0141    
0142 %    keyboard
0143    switch dirAxis
0144        case 'X'
0145            location= mean(xposition_elec(:,elecsUsed),2);
0146            depth= -abs(xposition_elec(:,elecsUsed(1))-xposition_elec(:,elecsUsed(2)))/depthRatio;
0147            electrodesLocation= unique(round(xposition_elec/depthPrecision)*depthPrecision);
0148        case 'Y'
0149            location= mean(yposition_elec(:,elecsUsed),2);
0150            depth= -abs(yposition_elec(:,elecsUsed(1))-yposition_elec(:,elecsUsed(2)))/depthRatio;
0151            electrodesLocation= unique(round(yposition_elec/depthPrecision)*depthPrecision);
0152        case 'Z'
0153            location= mean(zposition_elec(:,elecsUsed),2);
0154            depth= -abs(zposition_elec(:,elecsUsed(1))-zposition_elec(:,elecsUsed(2)))/depthRatio;
0155            electrodesLocation= unique(round(zposition_elec/depthPrecision)*depthPrecision);
0156        case 'R'
0157            location= mean(rposition_elec(:,elecsUsed),2);
0158            depth= -abs(rposition_elec(:,elecsUsed(1))-rposition_elec(:,elecsUsed(2)))/depthRatio;
0159            electrodesLocation= unique(round(rposition_elec/depthPrecision)*depthPrecision);
0160    end
0161     
0162    P= depth+1i*location;
0163    [Pu,iu,ju]= unique(round(P/depthPrecision)*depthPrecision);
0164    
0165    if length(Pu) < length(P)
0166        lu= location(iu);
0167        zu= depth(iu);
0168        du= iu*0;
0169        for i= 1:length(Pu)
0170            du(i)= mean(data(ju==i));
0171        end
0172        if length(fmdl.misc.sizepoint)>length(Pu)
0173            su= iu*0;
0174 %            misu= iu*0;
0175 %            masu= iu*0;
0176            for i= 1:length(Pu)
0177            su(i)= min(fmdl.misc.sizepoint(ju==i));
0178 %            misu(i)= min(fmdl.misc.sizepoint(ju==i));
0179 %            masu(i)= max(fmdl.misc.sizepoint(ju==i));
0180            end
0181            fmdl.misc.sizepoint= su;
0182 %            UU= [su misu masu masu-misu];
0183        end
0184    else
0185        lu= location;
0186        zu= depth;
0187        du= data;
0188    end
0189        
0190    scatter(lu,zu,fmdl.misc.sizepoint,(du),'filled','MarkerEdgeColor','k');
0191 %    hold on; plot(electrodesLocation,electrodesLocation*0,'x','Color',[0 0.5 0])
0192    xlabel('Distance (m)','fontsize',fs,'fontname','Times');
0193    ylabel('Pseudo-depth (m)','fontsize',fs,'fontname','Times')
0194 %     axis equal; axis tight;
0195    set(gca,'fontsize',fs,'fontname','Times')
0196 end
0197 
0198 
0199 function [depth,location,L,d]= plotPseudoSectionProfileSpecificUp(fmdl,data)
0200    fs= 20;
0201    [elec_posn,elecNumber]= electrodesPosition(fmdl);
0202    xposition_elec= reshape(elec_posn(elecNumber,1),[],4); 
0203    yposition_elec= reshape(elec_posn(elecNumber,2),[],4); 
0204    zposition_elec= reshape(elec_posn(elecNumber,3),[],4); zposition_elec= zposition_elec-min(zposition_elec(:));
0205    if isfield(fmdl.show_pseudosection,'minx') && isfield(fmdl.show_pseudosection,'maxx') && ...
0206        isfield(fmdl.show_pseudosection,'miny') && isfield(fmdl.show_pseudosection,'maxy')
0207        minx= fmdl.show_pseudosection.minx;
0208        maxx= fmdl.show_pseudosection.maxx;
0209        miny= fmdl.show_pseudosection.miny;
0210        maxy= fmdl.show_pseudosection.maxy;
0211    else
0212        xposition_elec= xposition_elec-min(xposition_elec(:));
0213        yposition_elec= yposition_elec-min(yposition_elec(:));
0214        minx= min(xposition_elec(:));
0215        maxx= max(xposition_elec(:));
0216        miny= min(yposition_elec(:));
0217        maxy= max(yposition_elec(:));
0218    end
0219    
0220    if isfield(fmdl.show_pseudosection,'xToDeduce') && isfield(fmdl.show_pseudosection,'yToDeduce')
0221        if strcmp(fmdl.show_pseudosection.xToDeduce,'minx')
0222            xToDeduce= minx;
0223        elseif strcmp(fmdl.show_pseudosection.xToDeduce,'maxx')
0224            xToDeduce= maxx;
0225        end
0226        if strcmp(fmdl.show_pseudosection.yToDeduce,'miny')
0227            yToDeduce= miny;
0228        elseif strcmp(fmdl.show_pseudosection.yToDeduce,'maxy')
0229            yToDeduce= maxy;
0230        end
0231    else
0232        xToDeduce= minx;
0233        yToDeduce= miny;
0234    end
0235    rposition_elec= sqrt((xposition_elec-xToDeduce).^2 + ...
0236        (yposition_elec-yToDeduce).^2);
0237      
0238 %    rposition_elec= sqrt((xposition_elec-min(xposition_elec(:))).^2 + ...
0239 %        (yposition_elec-max(yposition_elec(:))).^2);
0240    
0241    if isfield(fmdl.show_pseudosection,'elecsUsed')
0242        elecsUsed= fmdl.show_pseudosection.elecsUsed;
0243    else
0244        elecsUsed= 3:4;
0245    end
0246    
0247    if isfield(fmdl.show_pseudosection,'dirAxis')
0248        dirAxis= fmdl.show_pseudosection.dirAxis;
0249    else
0250        dirAxis= 'X';
0251    end
0252    
0253    if isfield(fmdl.show_pseudosection,'depthRatio')
0254        depthRatio= fmdl.show_pseudosection.depthRatio;
0255    else
0256        depthRatio= 3;
0257    end
0258    
0259     if isfield(fmdl.show_pseudosection,'depthPrecision')
0260        depthPrecision= fmdl.show_pseudosection.depthPrecision;
0261    else
0262        depthPrecision= 1;
0263     end
0264     
0265     if isfield(fmdl.show_pseudosection,'depthLevel')
0266        depthLevel= fmdl.show_pseudosection.depthLevel;
0267    else
0268        depthLevel= 0;
0269     end
0270    
0271    switch dirAxis
0272        case 'X'
0273            location= mean(xposition_elec(:,elecsUsed),2);
0274            depth= abs(xposition_elec(:,elecsUsed(1))-xposition_elec(:,elecsUsed(2)))/depthRatio;
0275            electrodesLocation= unique(round(xposition_elec/depthPrecision)*depthPrecision);
0276        case 'Y'
0277            location= mean(yposition_elec(:,elecsUsed),2);
0278            depth= abs(yposition_elec(:,elecsUsed(1))-yposition_elec(:,elecsUsed(2)))/depthRatio;
0279            electrodesLocation= unique(round(yposition_elec/depthPrecision)*depthPrecision);
0280        case 'Z'
0281            location= mean(zposition_elec(:,elecsUsed),2);
0282            depth= abs(zposition_elec(:,elecsUsed(1))-zposition_elec(:,elecsUsed(2)))/depthRatio;
0283            electrodesLocation= unique(round(zposition_elec/depthPrecision)*depthPrecision);
0284        case 'R'
0285            location= mean(rposition_elec(:,elecsUsed),2);%keyboard
0286            depth= abs(rposition_elec(:,elecsUsed(1))-rposition_elec(:,elecsUsed(2)))/depthRatio;
0287            electrodesLocation= unique(round(rposition_elec/depthPrecision)*depthPrecision);
0288    end
0289    L= (rposition_elec(:,2)-rposition_elec(:,1))/2;
0290    d= (rposition_elec(:,4)-rposition_elec(:,3))/2;
0291    
0292    
0293    depth= depth+depthLevel; 
0294    P= depth+1i*location;
0295    [Pu,iu,ju]= unique(round(P/depthPrecision)*depthPrecision);
0296 %    keyboard
0297    
0298    if length(Pu) < length(P)
0299        lu= location(iu);
0300        zu= depth(iu);
0301        du= iu*0;
0302        for i= 1:length(Pu)
0303            du(i)= mean(data(ju==i));
0304        end
0305        if length(fmdl.misc.sizepoint)>length(Pu)
0306            su= iu*0;
0307 %            misu= iu*0;
0308 %            masu= iu*0;
0309            for i= 1:length(Pu)
0310            su(i)= min(fmdl.misc.sizepoint(ju==i));
0311 %            misu(i)= min(fmdl.misc.sizepoint(ju==i));
0312 %            masu(i)= max(fmdl.misc.sizepoint(ju==i));
0313            end
0314            fmdl.misc.sizepoint= su;
0315 %            UU= [su misu masu masu-misu];
0316        end
0317    else
0318        lu= location;
0319        zu= depth;
0320        du= data;
0321    end
0322    
0323    
0324 %    keyboard
0325    scatter(lu,zu,fmdl.misc.sizepoint,(du),'filled','MarkerEdgeColor','k');
0326 %    hold on; plot(electrodesLocation,electrodesLocation*0,'x','Color',[0 0.5 0])
0327    xlabel('Distance (m)','fontsize',fs,'fontname','Times');
0328    ylabel('Pseudo-depth (m)','fontsize',fs,'fontname','Times')
0329 %     axis equal; axis tight;
0330    set(gca,'fontsize',fs,'fontname','Times')
0331 end
0332 
0333 
0334 
0335 function [zps,xps]= plotPseudoSectionProfileDown(fmdl,data)
0336    fs= 20;
0337    
0338    [elec_posn,elecNumber]= electrodesPosition(fmdl);
0339    xposition_elec= reshape(elec_posn(elecNumber,1),[],4);
0340    xps= mean(xposition_elec,2);
0341    
0342    AB= abs(xposition_elec(:,2)-xposition_elec(:,1));
0343    MN= abs(xposition_elec(:,4)-xposition_elec(:,3));
0344    AN= abs(xposition_elec(:,4)-xposition_elec(:,1));
0345    BM= abs(xposition_elec(:,3)-xposition_elec(:,2));
0346    AM= abs(xposition_elec(:,3)-xposition_elec(:,1));
0347    BN= abs(xposition_elec(:,4)-xposition_elec(:,2));
0348    
0349    a= MN; %max([AB MN AN BM AM BN ],[],2);
0350    zps= -a/3;
0351    
0352 %    a= abs(elecNumber(:,1)-elecNumber(:,2));
0353 %    de= elec_posn(1,1)-elec_posn(2,1);
0354 %
0355 %    % Identiy reciprocal data(elecNumber(:,1)-elecNumber(:,2)) > abs(elecNumber(:,3)-elecNumber(:,4))
0356 %    R= find(abs(elecNumber(:,1)-elecNumber(:,2)) < abs(elecNumber(:,3)-elecNumber(:,4)));
0357 %
0358 %    if ~isempty(R)
0359 %        xps(R)= (elec_posn(elecNumber(R,3),4)+elec_posn(elecNumber(R,4),4))/2;
0360 %        a(R)= abs(elecNumber(R,3)-elecNumber(R,4));
0361 %    end
0362      
0363    P= zps+1i*xps;
0364    [Pu,iu,ju]= unique(P);
0365    
0366    if length(Pu) < length(P)
0367        xu= xps(iu);
0368        zu= zps(iu);
0369        du= iu*0;
0370        for i= 1:length(Pu)
0371            du(i)= mean(data(ju==i));
0372        end
0373    else
0374        xu= xps;
0375        zu= zps;
0376        du= data;
0377    end
0378        
0379    scatter(xu,zu,fmdl.misc.sizepoint,(du),'filled','MarkerEdgeColor','k');
0380    xlabel('Distance (m)','fontsize',fs,'fontname','Times');
0381    ylabel('Pseudo-depth (m)','fontsize',fs,'fontname','Times')
0382     axis equal; axis tight;
0383    set(gca,'fontsize',fs,'fontname','Times')
0384 end
0385 
0386 function [xps,zps]= plotPseudoSectionProfileVert(fmdl,data)
0387    fs= 20;
0388    
0389    [elec_posn,elecNumber]= electrodesPosition(fmdl);
0390   
0391    zposition_elec= reshape(elec_posn(elecNumber,3),[],4);
0392    zps= mean(zposition_elec,2);
0393    
0394    AB= abs(zposition_elec(:,2)-zposition_elec(:,1));
0395    MN= abs(zposition_elec(:,4)-zposition_elec(:,3));
0396    AN= abs(zposition_elec(:,4)-zposition_elec(:,1));
0397    BM= abs(zposition_elec(:,3)-zposition_elec(:,2));
0398    AM= abs(zposition_elec(:,3)-zposition_elec(:,1));
0399    BN= abs(zposition_elec(:,4)-zposition_elec(:,2));
0400    
0401    a= max([AB MN AN BM AM BN ],[],2);
0402    xps= a/3;
0403       
0404    P= zps+1i*xps;
0405    [Pu,iu,ju]= unique(P);
0406    
0407    xu= xps(iu);
0408    zu= zps(iu);
0409    du= iu*0;
0410    for i= 1:length(Pu)
0411        du(i)= mean(data(ju==i));
0412    end
0413    
0414    scatter(xu,zu,fmdl.misc.sizepoint,(du),'filled','MarkerEdgeColor','k'); 
0415    xlabel('Pseudo distance (m)','fontsize',fs,'fontname','Times');
0416    if zps(1)<0
0417        ylabel('Depth (m)','fontsize',fs,'fontname','Times')
0418    else
0419        ylabel('Height (m)','fontsize',fs,'fontname','Times')
0420    end
0421    axis equal; axis tight;
0422    set(gca,'fontsize',fs,'fontname','Times')
0423 end
0424 
0425 function [r_point,th_point]= plotPseudoSectionCircularIn(fmdl,data)
0426 fs= 20;
0427 [elec_posn,elecNumber] = electrodesPosition(fmdl);
0428 [A,th_point,r,xc,yc] = polarPosition(elecNumber,elec_posn);
0429 a= max(A,[],2);
0430 r_point= a/2;
0431 r_point= (r-r_point/pi);%*9*pi*r/10+pi*r/10;
0432 
0433 [x_point,y_point]= pol2cart(th_point,r_point);
0434 x_point= x_point+xc; y_point= y_point+yc;
0435 
0436 P= x_point+1i*y_point;
0437 [Pu,iu,ju]= unique(P);
0438    
0439 xu= x_point(iu);
0440 zu= y_point(iu);
0441 du= iu*0;
0442 for i= 1:length(Pu)
0443     du(i)= mean(data(ju==i));
0444 end
0445    
0446 scatter(xu,zu,fmdl.misc.sizepoint,(du),'filled','MarkerEdgeColor','k');  
0447 xlabel('X (m)','fontsize',fs,'fontname','Times');
0448 ylabel('Y (m)','fontsize',fs,'fontname','Times')
0449 axis equal; axis tight;
0450 set(gca,'fontsize',fs,'fontname','Times')
0451 end
0452 
0453 function [dMN,th_point]= plotPseudoSectionCircularInVolcano(fmdl,data)
0454 fs= 20;
0455 [elec_posn,elecNumber] = electrodesPosition(fmdl);
0456 [xposition_elec,yposition_elec,zposition_elec] = electrodesPositionABMN(elecNumber,elec_posn);
0457 dMN= sqrt((xposition_elec(:,4)-xposition_elec(:,3)).^2 + ...
0458     (yposition_elec(:,4)-yposition_elec(:,3)).^2 + (zposition_elec(:,4)-zposition_elec(:,3)).^2);
0459 
0460 [A,r_bary,th_point,r,xc,yc,TH,R] = polarPosition(elecNumber,elec_posn);
0461 [th_bary,r_bary,th_bary34,r_bary34,r,xc,yc,barycenter_x34,barycenter_y34] = bary34(elecNumber,elec_posn);
0462 % a= A(:,2);
0463 % r_point= a;
0464 r= 400; 
0465 r_point= (r-r_bary)/(max(r-r_bary)-min(r-r_bary))*(r-5);%*9*pi*r/10+pi*r/10;
0466 r_point= r_point-min(r_point)+5;
0467 [x_point,y_point]= pol2cart(th_bary34,r_point);
0468 x_point= x_point+xc; y_point= y_point+yc;
0469 P= x_point+1i*y_point;
0470 [Pu,iu,ju]= unique(P);
0471 xu= x_point(iu); zu= y_point(iu); du= iu*0;
0472 for i= 1:length(Pu);     du(i)= mean(data(ju==i)); end
0473 figure; scatter(xu,zu,fmdl.misc.sizepoint,(du),'filled','MarkerEdgeColor','k');  
0474 text(xu(iu==1),zu(iu==1),'1')
0475 hold on; plot(xposition_elec,yposition_elec,'kp')
0476 xlabel('X (m)','fontsize',fs,'fontname','Times');
0477 ylabel('Y (m)','fontsize',fs,'fontname','Times')
0478 axis equal; axis tight;
0479 set(gca,'fontsize',fs,'fontname','Times')
0480 end
0481 
0482 function [r_point,th_point]= plotPseudoSectionCircularOut(fmdl,data)
0483 fs= 20;
0484 [elec_posn,elecNumber] = electrodesPosition(fmdl);
0485 [A,r_bary,th_point,r,xc,yc] = polarPosition(elecNumber,elec_posn);
0486 a= max(A,[],2);
0487 r_point= a/2;
0488 r_point= (r+r_point);
0489 [x_point,y_point]= pol2cart(th_point,r_point+r);
0490 x_point= x_point+xc; y_point= y_point+yc;
0491 
0492 P= x_point+1i*y_point;
0493 [Pu,iu,ju]= unique(P);
0494    
0495 xu= x_point(iu);
0496 zu= y_point(iu);
0497 du= iu*0;
0498 for i= 1:length(Pu)
0499     du(i)= mean(data(ju==i));
0500 end
0501    
0502 scatter(xu,zu,fmdl.misc.sizepoint,(du),'filled','MarkerEdgeColor','k');  colorbar
0503 xlabel('X (m)','fontsize',fs,'fontname','Times');
0504 ylabel('Y (m)','fontsize',fs,'fontname','Times')
0505 axis equal; axis tight;
0506 set(gca,'fontsize',fs,'fontname','Times')
0507 end
0508 
0509 
0510 function [elec_posn,elecNumber] = electrodesPosition(fmdl)
0511     stim= fmdl.stimulation;
0512     stimulationMatrix= [];
0513     for i = 1:length(stim);
0514         nmp= size(stim(i).meas_pattern, 1);
0515         [idxIN,idxJN]= find(stim(i).meas_pattern<0);
0516         [idxIP,idxJP]= find(stim(i).meas_pattern>0);
0517         stimulationMatrix= [stimulationMatrix; [ 0*(1:nmp)'+find(stim(i).stim_pattern<0) ...
0518             0*(1:nmp)'+find(stim(i).stim_pattern>0)] idxJN(idxIN) idxJP(idxIP)];
0519     end
0520 %    A=  stimulationMatrix(:,1);
0521 %    A=  stimulationMatrix(:,1);
0522 %    A=  stimulationMatrix(:,1);
0523 %    A=  stimulationMatrix(:,1);
0524 %
0525 %    A= zeros(length(fmdl.stimulation),1);
0526 %    B= zeros(length(fmdl.stimulation),1);
0527 %    M= zeros(length(fmdl.stimulation),1);
0528 %    N= zeros(length(fmdl.stimulation),1);
0529 %    for i=1:length(fmdl.stimulation)
0530 %        A(i)= find(fmdl.stimulation(1,i).stim_pattern<0);
0531 %        B(i)= find(fmdl.stimulation(1,i).stim_pattern>0);
0532 %        M(i)= find(fmdl.stimulation(1,i).meas_pattern<0);
0533 %        N(i)= find(fmdl.stimulation(1,i).meas_pattern>0);
0534 %    end
0535 %    elecNumber= [A B M N];
0536    
0537    elecNumber= stimulationMatrix;
0538    
0539    elec_posn= zeros(length(fmdl.electrode),size(fmdl.nodes,2));
0540    for i=1:length(fmdl.electrode)
0541    elec_posn(i,:)= mean(fmdl.nodes(fmdl.electrode(1,i).nodes,:),1);
0542    end
0543    
0544 end
0545 
0546 function [xposition_elec,yposition_elec,zposition_elec] = electrodesPositionABMN(elecNumber,elec_posn)
0547 xposition_elec= reshape(elec_posn(elecNumber,1),[],4);
0548 yposition_elec= reshape(elec_posn(elecNumber,2),[],4);
0549 zposition_elec= reshape(elec_posn(elecNumber,3),[],4);
0550 end
0551 
0552 function [A,r_bary,th_bary,r,xc,yc,TH,R] = polarPosition(elecNumber,elec_posn)
0553 TH= elecNumber*0;
0554 R= elecNumber*0;
0555 
0556 xposition_elec= reshape(elec_posn(elecNumber,1),[],4);
0557 yposition_elec= reshape(elec_posn(elecNumber,2),[],4);
0558 rx= (max(xposition_elec(:))-min(xposition_elec(:)))/2;
0559 ry= (max(yposition_elec(:))-min(yposition_elec(:)))/2;
0560 rmean= (rx+ry)/2;
0561 
0562 xc= mean(elec_posn(:,1)); yc= mean(elec_posn(:,2)); r = rmean;
0563 
0564 barycenter_x= mean(xposition_elec,2);
0565 barycenter_y= mean(yposition_elec,2);
0566 
0567 [th_bary,r_bary]= cart2pol(barycenter_x-xc,barycenter_y-yc);
0568 
0569 for i= 1:4
0570     [TH(:,i),R(:,i)]= cart2pol(elec_posn(elecNumber(:,i),1)-xc,elec_posn(elecNumber(:,i),2)-yc);
0571 end 
0572 
0573 Arc_length_AB= (TH(:,2)-TH(:,1));
0574 idx_discontinuity= find(TH(:,2)>TH(:,1) & TH(:,2)<0);
0575 Arc_length_AB(idx_discontinuity)= (-TH(idx_discontinuity,1)+TH(idx_discontinuity,2));
0576 idx_discontinuity= find(TH(:,1)>TH(:,2));
0577 Arc_length_AB(idx_discontinuity)= (TH(idx_discontinuity,1)-TH(idx_discontinuity,2));
0578 idx_discontinuity= find(TH(:,1)>TH(:,2) & TH(:,1)<0);
0579 Arc_length_AB(idx_discontinuity)= (-TH(idx_discontinuity,2)+TH(idx_discontinuity,1));
0580 Arc_length_AB(Arc_length_AB>=pi+0.005)= 2*pi-Arc_length_AB(Arc_length_AB>=pi+0.005);
0581 Arc_length_AB= Arc_length_AB.*(R(:,1)+R(:,2))/2;
0582 
0583 Arc_length_MN= (TH(:,4)-TH(:,3));
0584 idx_discontinuity= find((TH(:,4)>TH(:,3) & TH(:,4)<0));
0585 Arc_length_MN(idx_discontinuity)= -TH(idx_discontinuity,3)+TH(idx_discontinuity,4);
0586 idx_discontinuity= find((TH(:,4)<TH(:,3) & TH(:,4)<0));
0587 Arc_length_MN(idx_discontinuity)= +TH(idx_discontinuity,3)-TH(idx_discontinuity,4);
0588 idx_discontinuity= find((TH(:,3)>TH(:,4)) & TH(:,3)>0);
0589 Arc_length_MN(idx_discontinuity)= TH(idx_discontinuity,3)-TH(idx_discontinuity,4);
0590 % idx_discontinuity= find((TH(:,3)>TH(:,4)) &  TH(:,3)<0);
0591 % Arc_length_MN(idx_discontinuity)= -TH(idx_discontinuity,3)+TH(idx_discontinuity,4);
0592 Arc_length_MN(Arc_length_MN>=pi+0.005)= 2*pi-Arc_length_MN(Arc_length_MN>=pi+0.005);
0593 Arc_length_MN= Arc_length_MN.*(R(:,3)+R(:,4))/2;
0594 
0595 Arc_length_AM= (TH(:,3)-TH(:,1));
0596 idx_discontinuity= find((TH(:,3)>TH(:,1) & TH(:,3)<0));
0597 Arc_length_AM(idx_discontinuity)= -TH(idx_discontinuity,1)+TH(idx_discontinuity,3);
0598 idx_discontinuity= find((TH(:,1)> TH(:,3)));
0599 Arc_length_AM(idx_discontinuity)= TH(idx_discontinuity,1)-TH(idx_discontinuity,3);
0600 idx_discontinuity= find((TH(:,1)> TH(:,3) & TH(:,1) <0));
0601 Arc_length_AM(idx_discontinuity)= -TH(idx_discontinuity,1)+TH(idx_discontinuity,3);
0602 Arc_length_AM(Arc_length_AM>=pi+0.005)= 2*pi-Arc_length_AM(Arc_length_AM>=pi+0.005);
0603 Arc_length_AM= Arc_length_AM.*(R(:,1)+R(:,3))/2;
0604 
0605 Arc_length_AN= (TH(:,4)-TH(:,1));
0606 idx_discontinuity= find((TH(:,4)>TH(:,1)<0 & TH(:,4)<0));
0607 Arc_length_AN(idx_discontinuity)= -TH(idx_discontinuity,1)+TH(idx_discontinuity,4);
0608 idx_discontinuity= find((TH(:,1)>TH(:,4)));
0609 Arc_length_AN(idx_discontinuity)= TH(idx_discontinuity,1)-TH(idx_discontinuity,4);
0610 idx_discontinuity= find((TH(:,1)>TH(:,4) & TH(:,1)<0));
0611 Arc_length_AN(idx_discontinuity)= TH(idx_discontinuity,4)-TH(idx_discontinuity,1);
0612 Arc_length_AN(Arc_length_AN>=pi+0.005)= 2*pi-Arc_length_AN(Arc_length_AN>=pi+0.005);
0613 Arc_length_AN= Arc_length_AN.*(R(:,1)+R(:,4))/2;
0614 
0615 Arc_length_NB= (TH(:,4)-TH(:,2));
0616 idx_discontinuity= find((TH(:,4)>TH(:,2)));
0617 Arc_length_NB(idx_discontinuity)= -TH(idx_discontinuity,2)+TH(idx_discontinuity,4);
0618 idx_discontinuity= find((TH(:,2)> TH(:,4)));
0619 Arc_length_NB(idx_discontinuity)= -TH(idx_discontinuity,4)+TH(idx_discontinuity,2);
0620 idx_discontinuity= find((TH(:,2)> TH(:,4) &  TH(:,4)<0));
0621 Arc_length_NB(idx_discontinuity)= -TH(idx_discontinuity,4)+TH(idx_discontinuity,2);
0622 Arc_length_NB(Arc_length_NB>=pi+0.005)= 2*pi-Arc_length_NB(Arc_length_NB>=pi+0.005);
0623 Arc_length_NB= Arc_length_NB.*(R(:,2)+R(:,4))/2;
0624 
0625 Arc_length_BM= (TH(:,3)-TH(:,2));
0626 idx_discontinuity= find((TH(:,3)>TH(:,2) & TH(:,2)<0));
0627 Arc_length_BM(idx_discontinuity)= -TH(idx_discontinuity,2)+TH(idx_discontinuity,3);
0628 idx_discontinuity= find((TH(:,2)> TH(:,3)));
0629 Arc_length_BM(idx_discontinuity)= -TH(idx_discontinuity,3)+TH(idx_discontinuity,2);
0630 idx_discontinuity= find((TH(:,2)> TH(:,3)& TH(:,3)<0));
0631 Arc_length_BM(idx_discontinuity)= TH(idx_discontinuity,3)-TH(idx_discontinuity,2);
0632 Arc_length_BM(Arc_length_BM>=pi+0.005)= 2*pi-Arc_length_BM(Arc_length_BM>=pi+0.005);
0633 Arc_length_BM= Arc_length_BM.*(R(:,2)+R(:,3))/2;
0634 
0635 A= [Arc_length_AB Arc_length_MN Arc_length_AM Arc_length_AN Arc_length_NB Arc_length_BM];
0636 end
0637 
0638 
0639 
0640 function [th_bary,r_bary,th_bary34,r_bary34,r,xc,yc,barycenter_x34,barycenter_y34] = bary34(elecNumber,elec_posn)
0641 xposition_elec= reshape(elec_posn(elecNumber,1),[],4);
0642 yposition_elec= reshape(elec_posn(elecNumber,2),[],4);
0643 rx= (max(xposition_elec(:))-min(xposition_elec(:)))/2;
0644 ry= (max(yposition_elec(:))-min(yposition_elec(:)))/2;
0645 rmean= (rx+ry)/2;
0646 
0647 xc= mean(elec_posn(:,1)); yc= mean(elec_posn(:,2)); r = rmean;
0648 
0649 barycenter_x34= mean(xposition_elec(:,3:4),2);
0650 barycenter_y34= mean(yposition_elec(:,3:4),2);
0651 [th_bary34,r_bary34]= cart2pol(barycenter_x34-xc,barycenter_y34-yc);
0652 
0653 barycenter_x= mean(xposition_elec,2);
0654 barycenter_y= mean(yposition_elec,2);
0655 [th_bary,r_bary]= cart2pol(barycenter_x-xc,barycenter_y-yc);
0656 
0657 
0658 end
0659 
0660 function do_unit_test
0661 shape_str = ['solid top    = plane(0,0,0;0,1,0);\n' ...
0662              'solid mainobj= top and orthobrick(-100,-200,-100;410,10,100) -maxh=20.0;\n'];
0663 e0 = linspace(0,310,64)';
0664 elec_pos = [e0,0*e0,0*e0,1+0*e0,0*e0,0*e0];
0665 elec_shape= [0.1,0.1,1];
0666 elec_obj = 'top';
0667 fmdl = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0668 
0669 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];
0670 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];
0671 fmdl.stimulation= stim_pattern_geophys( 64, 'DipoleDipole', {'spacings', spacing,'multiples',multiples} );
0672 
0673 img1= mk_image(fmdl,1);
0674 vh1= fwd_solve(img1);
0675 normalisation= 1./vh1.meas;
0676 I= speye(length(normalisation));
0677 I(1:size(I,1)+1:size(I,1)*size(I,1))= normalisation;
0678 
0679 img = mk_image(fmdl,0+ mk_c2f_circ_mapping(fmdl,[100;-50;0;50])*100);
0680 img.elem_data(img.elem_data==0)= 0.1;
0681 dd  = fwd_solve(img);
0682 subplot(221);
0683 show_pseudosection( fmdl, I*dd.meas )
0684 
0685 subplot(222);
0686 fmdl.show_pseudosection.orientation = 'horizontaldownward';
0687 show_pseudosection( fmdl, I*dd.meas )
0688 
0689 subplot(223);
0690 fmdl.show_pseudosection.orientation = 'vertical';
0691 show_pseudosection( fmdl, I*dd.meas )
0692 
0693 end

Generated on Tue 12-May-2015 16:00:42 by m2html © 2005