show_current

 SHOW_CURRENT: show current or other quantity defined
  on nodes onto the image

 show_current( img, volt )
   img -> img object 
   volt-> voltage on nodes (if not specified, img is solved
      via fwd_solve, or the value on node_data is used)

 Without specified values, show_current will
   create one current vector for each element
 The points are specified as either
   img.fwd_model.mdl_slice_mapper.npx   - number of points in horizontal direction
   img.fwd_model.mdl_slice_mapper.npy   - number of points in vertical 
    or
   img.fwd_model.mdl_slice_mapper.x_pts - vector of points in horizontal direction
   img.fwd_model.mdl_slice_mapper.y_pts - vector of points in vertical
    or
   img.fwd_model.mdl_slice_mapper.npoints    - number of points across the larger
                                               dimension of the model (square)
    or
   img.fwd_model.mdl_slice_mapper.resolution - number of points per unit

 For 3D models, the slice may be specified as (see mdl_slice_mapper for information)
   img.fwd_model.mdl_slice_mapper.level 
 
 If an output is specified, then no image is draws
 
 Examples:
   img = mk_image( mk_common_model('b2c2',8));
   show_current(img); 
 OR
   img = mk_image( mk_common_model('b2c2',8));
   img.fwd_model.mdl_slice_mapper.npx  = 64;
   img.fwd_model.mdl_slice_mapper.npy  = 32;
   show_current(img); 
 OR
   img = mk_image( mk_common_model('b2c2',8));
   q = show_current(img); 
   quiver(q.xp,q.yp, q.xc,q.yc, 2,'b','LineWidth',2);

0001 function quiv = show_current( img, vv )
0002 % SHOW_CURRENT: show current or other quantity defined
0003 %  on nodes onto the image
0004 %
0005 % show_current( img, volt )
0006 %   img -> img object
0007 %   volt-> voltage on nodes (if not specified, img is solved
0008 %      via fwd_solve, or the value on node_data is used)
0009 %
0010 % Without specified values, show_current will
0011 %   create one current vector for each element
0012 % The points are specified as either
0013 %   img.fwd_model.mdl_slice_mapper.npx   - number of points in horizontal direction
0014 %   img.fwd_model.mdl_slice_mapper.npy   - number of points in vertical
0015 %    or
0016 %   img.fwd_model.mdl_slice_mapper.x_pts - vector of points in horizontal direction
0017 %   img.fwd_model.mdl_slice_mapper.y_pts - vector of points in vertical
0018 %    or
0019 %   img.fwd_model.mdl_slice_mapper.npoints    - number of points across the larger
0020 %                                               dimension of the model (square)
0021 %    or
0022 %   img.fwd_model.mdl_slice_mapper.resolution - number of points per unit
0023 %
0024 % For 3D models, the slice may be specified as (see mdl_slice_mapper for information)
0025 %   img.fwd_model.mdl_slice_mapper.level
0026 %
0027 % If an output is specified, then no image is draws
0028 %
0029 % Examples:
0030 %   img = mk_image( mk_common_model('b2c2',8));
0031 %   show_current(img);
0032 % OR
0033 %   img = mk_image( mk_common_model('b2c2',8));
0034 %   img.fwd_model.mdl_slice_mapper.npx  = 64;
0035 %   img.fwd_model.mdl_slice_mapper.npy  = 32;
0036 %   show_current(img);
0037 % OR
0038 %   img = mk_image( mk_common_model('b2c2',8));
0039 %   q = show_current(img);
0040 %   quiver(q.xp,q.yp, q.xc,q.yc, 2,'b','LineWidth',2);
0041 
0042 
0043 % (C) 2010 Andy Adler. License: GPL version 2 or version 3
0044 % $Id: show_current.m 6340 2022-04-21 21:12:36Z bgrychtol $
0045 
0046 if ischar(img) && strcmp(img,'UNIT_TEST'); do_unit_test; return; end
0047 
0048 
0049 if nargin==1;
0050   e_curr = calc_elem_current( img );
0051 else 
0052   e_curr = calc_elem_current( img, vv);
0053 end
0054 
0055 dims = size(img.fwd_model.nodes,2);
0056 
0057 elemcur = [zeros(1,dims); e_curr ];
0058 try
0059    if strcmp(img.calc_colours.component, 'real')
0060        elemcur = real(elemcur);
0061    end
0062    if strcmp(img.calc_colours.component, 'imag')
0063        elemcur = imag(elemcur);
0064    end
0065 end  % if no calc_colours.component field, do nothing
0066 
0067 if isfield(img.fwd_model, 'mdl_slice_mapper');
0068    if dims == 2;
0069       img.fwd_model.mdl_slice_mapper.level = [inf,inf,0];
0070    end
0071    elem_ptr = mdl_slice_mapper( img.fwd_model, 'elem' );
0072    szep = size(elem_ptr);
0073 
0074    %[xp,yp] = grid_the_space( img.fwd_model);
0075    xy = mdl_slice_mapper( img.fwd_model, 'get_points');
0076    xp = xy{1}; yp = xy{2};
0077 
0078    xc = reshape( elemcur(elem_ptr+1,1), szep);
0079    yc = reshape( elemcur(elem_ptr+1,2), szep);
0080    if dims==3
0081       zc = reshape( elemcur(elem_ptr+1,3), szep);
0082    end
0083 else 
0084    pp = interp_mesh( img.fwd_model);
0085    xp = pp(:,1);
0086    yp= pp(:,2);
0087    xc = elemcur(2:end,1);
0088    yc = elemcur(2:end,2);
0089    if dims==3
0090       zc = elemcur(2:end,3);
0091    end
0092 end
0093 
0094 quiv.xp = xp; 
0095 quiv.yp = yp; 
0096 quiv.xc = xc; 
0097 quiv.yc = yc; 
0098 if dims==3
0099    quiv.zc = zc; 
0100 end
0101 if nargout == 0
0102    quiver(xp,yp,xc,yc,2,'k');
0103 end
0104 
0105 function vv = fix_dim(vv)
0106     if size(vv,1) == 1
0107         vv = vv';
0108     end
0109 
0110 function  [x,y] = grid_the_space( fmdl);
0111 
0112   xspace = []; yspace = [];
0113   try 
0114      xspace = fmdl.mdl_slice_mapper.x_pts;
0115      yspace = fmdl.mdl_slice_mapper.y_pts;
0116   end
0117 
0118   if isempty(xspace)
0119      npx  = fmdl.mdl_slice_mapper.npx;
0120      npy  = fmdl.mdl_slice_mapper.npy;
0121 
0122      xmin = min(fmdl.nodes(:,1));    xmax = max(fmdl.nodes(:,1));
0123      xmean= mean([xmin,xmax]); xrange= xmax-xmin;
0124 
0125      ymin = min(fmdl.nodes(:,2));    ymax = max(fmdl.nodes(:,2));
0126      ymean= mean([ymin,ymax]); yrange= ymax-ymin;
0127 
0128      range= max([xrange, yrange]);
0129      xspace = linspace( xmean - range*0.5, xmean + range*0.5, npx );
0130      yspace = linspace( ymean + range*0.5, ymean - range*0.5, npy );
0131   end
0132   if size(xspace,2) == 1
0133       [x,y]=meshgrid( xspace, yspace );
0134   else
0135       x= xspace;
0136       y= yspace;
0137   end
0138 
0139 function do_unit_test
0140    fmdl.nodes = [0,0;0,1;1,0;1,1];
0141    fmdl.elems = [1,2,3;2,3,4];
0142    fmdl.electrode(1).nodes = [1,2]; fmdl.electrode(1).z_contact = 0.01;
0143    fmdl.electrode(2).nodes = [3,4]; fmdl.electrode(2).z_contact = 0.01;
0144    fmdl.gnd_node = 1;
0145    fmdl.stimulation(1).stim_pattern = [1;-1];
0146    fmdl.stimulation(1).meas_pattern = [1,-1];
0147    fmdl.solve = @fwd_solve_1st_order;
0148    fmdl.system_mat = @system_mat_1st_order;
0149    fmdl.type = 'fwd_model';
0150    fmdl.normalize_measurements= 0;
0151    img = mk_image(fmdl,[1,1]); 
0152    img.fwd_solve.get_all_meas = 1;
0153 
0154    img.fwd_model.mdl_slice_mapper.npx = 6;
0155    img.fwd_model.mdl_slice_mapper.npy = 6;
0156    show_current(img);
0157 
0158    imdl= mk_common_model('d2d2c',8);
0159    img = calc_jacobian_bkgnd( imdl );
0160    img.fwd_model.mdl_slice_mapper.npoints = 64;
0161    show_current(img);
0162 
0163    img.fwd_model = rmfield(img.fwd_model, 'mdl_slice_mapper');
0164    img.fwd_model.mdl_slice_mapper.x_pts = linspace(0,1,62);
0165    img.fwd_model.mdl_slice_mapper.y_pts = linspace(0,1,56);
0166    img.fwd_model.mdl_slice_mapper.level = [inf,inf,0];
0167 
0168    img.fwd_solve.get_all_meas = 1;
0169    vh = fwd_solve(img);
0170    show_current(img, vh.volt(:,1));