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

 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 %
0019 % For 3D models, the slice may be specified as (see mdl_slice_mapper for information)
0020 %   img.fwd_model.mdl_slice_mapper.level
0021 %
0022 % If an output is specified, then no image is draws
0023 %
0024 % Examples:
0025 %   img = mk_image( mk_common_model('b2c2',8));
0026 %   show_current(img);
0027 % OR
0028 %   img = mk_image( mk_common_model('b2c2',8));
0029 %   img.fwd_model.mdl_slice_mapper.npx  = 64;
0030 %   img.fwd_model.mdl_slice_mapper.npy  = 32;
0031 %   show_current(img);
0032 % OR
0033 %   img = mk_image( mk_common_model('b2c2',8));
0034 %   q = show_current(img);
0035 %   quiver(q.xp,q.yp, q.xc,q.yc, 2,'b','LineWidth',2);
0036 
0037 
0038 % (C) 2010 Andy Adler. License: GPL version 2 or version 3
0039 % $Id: show_current.m 5541 2017-06-15 11:34:31Z aadler $
0040 
0041 if ischar(img) && strcmp(img,'UNIT_TEST'); do_unit_test; return; end
0042 
0043 
0044 if nargin==1;
0045   e_curr = calc_elem_current( img );
0046 else 
0047   e_curr = calc_elem_current( img, vv);
0048 end
0049 
0050 dims = size(img.fwd_model.nodes,2);
0051 
0052 elemcur = [zeros(1,dims); e_curr ];
0053 try
0054    if strcmp(img.calc_colours.component, 'real')
0055        elemcur = real(elemcur);
0056    end
0057    if strcmp(img.calc_colours.component, 'imag')
0058        elemcur = imag(elemcur);
0059    end
0060 end  % if no calc_colours.component field, do nothing
0061 
0062 if isfield(img.fwd_model, 'mdl_slice_mapper');
0063    if dims == 2;
0064       img.fwd_model.mdl_slice_mapper.level = [inf,inf,0];
0065    end
0066    elem_ptr = mdl_slice_mapper( img.fwd_model, 'elem' );
0067    szep = size(elem_ptr);
0068 
0069    [xp,yp] = grid_the_space( img.fwd_model);
0070 
0071    xc = reshape( elemcur(elem_ptr+1,1), szep);
0072    yc = reshape( elemcur(elem_ptr+1,2), szep);
0073    if dims==3
0074       zc = reshape( elemcur(elem_ptr+1,3), szep);
0075    end
0076 else 
0077    pp = interp_mesh( img.fwd_model);
0078    xp = pp(:,1);
0079    yp= pp(:,2);
0080    xc = elemcur(2:end,1);
0081    yc = elemcur(2:end,2);
0082    if dims==3
0083       zc = elemcur(2:end,3);
0084    end
0085 end
0086 
0087 quiv.xp = xp; 
0088 quiv.yp = yp; 
0089 quiv.xc = xc; 
0090 quiv.yc = yc; 
0091 if dims==3
0092    quiv.zc = zc; 
0093 end
0094 if nargout == 0
0095    quiver(xp,yp,xc,yc,2,'k');
0096 end
0097 
0098 function vv = fix_dim(vv)
0099     if size(vv,1) == 1
0100         vv = vv';
0101     end
0102 
0103 function  [x,y] = grid_the_space( fmdl);
0104 
0105   xspace = []; yspace = [];
0106   try 
0107      xspace = fmdl.mdl_slice_mapper.x_pts;
0108      yspace = fmdl.mdl_slice_mapper.y_pts;
0109   end
0110 
0111   if isempty(xspace)
0112      npx  = fmdl.mdl_slice_mapper.npx;
0113      npy  = fmdl.mdl_slice_mapper.npy;
0114 
0115      xmin = min(fmdl.nodes(:,1));    xmax = max(fmdl.nodes(:,1));
0116      xmean= mean([xmin,xmax]); xrange= xmax-xmin;
0117 
0118      ymin = min(fmdl.nodes(:,2));    ymax = max(fmdl.nodes(:,2));
0119      ymean= mean([ymin,ymax]); yrange= ymax-ymin;
0120 
0121      range= max([xrange, yrange]);
0122      xspace = linspace( xmean - range*0.5, xmean + range*0.5, npx );
0123      yspace = linspace( ymean + range*0.5, ymean - range*0.5, npy );
0124   end
0125   if size(xspace,2) == 1
0126       [x,y]=meshgrid( xspace, yspace );
0127   else
0128       x= xspace;
0129       y= yspace;
0130   end
0131 
0132 function do_unit_test
0133    fmdl.nodes = [0,0;0,1;1,0;1,1];
0134    fmdl.elems = [1,2,3;2,3,4];
0135    fmdl.electrode(1).nodes = [1,2]; fmdl.electrode(1).z_contact = 0.01;
0136    fmdl.electrode(2).nodes = [3,4]; fmdl.electrode(2).z_contact = 0.01;
0137    fmdl.gnd_node = 1;
0138    fmdl.stimulation(1).stim_pattern = [1;-1];
0139    fmdl.stimulation(1).meas_pattern = [1,-1];
0140    fmdl.solve = @fwd_solve_1st_order;
0141    fmdl.system_mat = @system_mat_1st_order;
0142    fmdl.type = 'fwd_model';
0143    fmdl.normalize_measurements= 0;
0144    img = mk_image(fmdl,[1,1]); 
0145    img.fwd_solve.get_all_meas = 1;
0146 
0147    img.fwd_model.mdl_slice_mapper.npx = 6;
0148    img.fwd_model.mdl_slice_mapper.npy = 6;
0149    show_current(img);
0150 
0151    imdl= mk_common_model('d2d2c',8);
0152    img = calc_jacobian_bkgnd( imdl );
0153    img.fwd_model.mdl_slice_mapper.npx = 64;
0154    img.fwd_model.mdl_slice_mapper.npy = 64;
0155    show_current(img);
0156 
0157    img.fwd_model.mdl_slice_mapper.x_pts = linspace(0,1,62);
0158    img.fwd_model.mdl_slice_mapper.y_pts = linspace(0,1,56);
0159    img.fwd_model.mdl_slice_mapper.level = [inf,inf,0];
0160 
0161    img.fwd_solve.get_all_meas = 1;
0162    vh = fwd_solve(img);
0163    show_current(img, vh.volt(:,1));