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);

 See also:
  calc_elem_current

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 % See also:
0043 %  calc_elem_current
0044 
0045 
0046 % (C) 2010 Andy Adler. License: GPL version 2 or version 3
0047 % $Id: show_current.m 6941 2024-06-18 11:58:10Z aadler $
0048 
0049 if ischar(img) && strcmp(img,'UNIT_TEST'); do_unit_test; return; end
0050 
0051 
0052 if nargin==1;
0053   e_curr = calc_elem_current( img );
0054 else 
0055   e_curr = calc_elem_current( img, vv);
0056 end
0057 
0058 dims = size(img.fwd_model.nodes,2);
0059 
0060 elemcur = [zeros(1,dims); e_curr ];
0061 try
0062    if strcmp(img.calc_colours.component, 'real')
0063        elemcur = real(elemcur);
0064    end
0065    if strcmp(img.calc_colours.component, 'imag')
0066        elemcur = imag(elemcur);
0067    end
0068 end  % if no calc_colours.component field, do nothing
0069 
0070 if isfield(img.fwd_model, 'mdl_slice_mapper');
0071    if dims == 2;
0072       img.fwd_model.mdl_slice_mapper.level = [inf,inf,0];
0073    end
0074    elem_ptr = mdl_slice_mapper( img.fwd_model, 'elem' );
0075    szep = size(elem_ptr);
0076 
0077    %[xp,yp] = grid_the_space( img.fwd_model);
0078    xy = mdl_slice_mapper( img.fwd_model, 'get_points');
0079    xp = xy{1}; yp = xy{2};
0080 
0081    xc = reshape( elemcur(elem_ptr+1,1), szep);
0082    yc = reshape( elemcur(elem_ptr+1,2), szep);
0083    if dims==3
0084       zc = reshape( elemcur(elem_ptr+1,3), szep);
0085    end
0086 else 
0087    pp = interp_mesh( img.fwd_model);
0088    xp = pp(:,1);
0089    yp= pp(:,2);
0090    xc = elemcur(2:end,1);
0091    yc = elemcur(2:end,2);
0092    if dims==3
0093       zc = elemcur(2:end,3);
0094    end
0095 end
0096 
0097 quiv.xp = xp; 
0098 quiv.yp = yp; 
0099 quiv.xc = xc; 
0100 quiv.yc = yc; 
0101 if dims==3
0102    quiv.zc = zc; 
0103 end
0104 if nargout == 0
0105    quiver(xp,yp,xc,yc,2,'k');
0106 end
0107 
0108 function vv = fix_dim(vv)
0109     if size(vv,1) == 1
0110         vv = vv';
0111     end
0112 
0113 function  [x,y] = grid_the_space( fmdl);
0114 
0115   xspace = []; yspace = [];
0116   try 
0117      xspace = fmdl.mdl_slice_mapper.x_pts;
0118      yspace = fmdl.mdl_slice_mapper.y_pts;
0119   end
0120 
0121   if isempty(xspace)
0122      npx  = fmdl.mdl_slice_mapper.npx;
0123      npy  = fmdl.mdl_slice_mapper.npy;
0124 
0125      xmin = min(fmdl.nodes(:,1));    xmax = max(fmdl.nodes(:,1));
0126      xmean= mean([xmin,xmax]); xrange= xmax-xmin;
0127 
0128      ymin = min(fmdl.nodes(:,2));    ymax = max(fmdl.nodes(:,2));
0129      ymean= mean([ymin,ymax]); yrange= ymax-ymin;
0130 
0131      range= max([xrange, yrange]);
0132      xspace = linspace( xmean - range*0.5, xmean + range*0.5, npx );
0133      yspace = linspace( ymean + range*0.5, ymean - range*0.5, npy );
0134   end
0135   if size(xspace,2) == 1
0136       [x,y]=meshgrid( xspace, yspace );
0137   else
0138       x= xspace;
0139       y= yspace;
0140   end
0141 
0142 function do_unit_test
0143    fmdl.nodes = [0,0;0,1;1,0;1,1];
0144    fmdl.elems = [1,2,3;2,3,4];
0145    fmdl.electrode(1).nodes = [1,2]; fmdl.electrode(1).z_contact = 0.01;
0146    fmdl.electrode(2).nodes = [3,4]; fmdl.electrode(2).z_contact = 0.01;
0147    fmdl.gnd_node = 1;
0148    fmdl.stimulation(1).stim_pattern = [1;-1];
0149    fmdl.stimulation(1).meas_pattern = [1,-1];
0150    fmdl.solve = @fwd_solve_1st_order;
0151    fmdl.system_mat = @system_mat_1st_order;
0152    fmdl.type = 'fwd_model';
0153    fmdl.normalize_measurements= 0;
0154    img = mk_image(fmdl,[1,1]); 
0155    img.fwd_solve.get_all_meas = 1;
0156 
0157    img.fwd_model.mdl_slice_mapper.npx = 6;
0158    img.fwd_model.mdl_slice_mapper.npy = 6;
0159    show_current(img);
0160 
0161    imdl= mk_common_model('d2d2c',8);
0162    img = calc_jacobian_bkgnd( imdl );
0163    img.fwd_model.mdl_slice_mapper.npoints = 64;
0164    show_current(img);
0165 
0166    img.fwd_model = rmfield(img.fwd_model, 'mdl_slice_mapper');
0167    img.fwd_model.mdl_slice_mapper.x_pts = linspace(0,1,62);
0168    img.fwd_model.mdl_slice_mapper.y_pts = linspace(0,1,56);
0169    img.fwd_model.mdl_slice_mapper.level = [inf,inf,0];
0170 
0171    img.fwd_solve.get_all_meas = 1;
0172    vh = fwd_solve(img);
0173    show_current(img, vh.volt(:,1));