% COMSOL Multiphysics Model M-file
% Generated by COMSOL 3.5a (COMSOL 3.5.0.603, $Date: 2008/12/03 17:02:19 $)
% Some geometry objects are stored in a separate file.
% The name of this file is given by the variable 'flbinaryfile'.

set (0, 'DefaultFigureWindowStyle', 'docked')
clc            % clear the command window
clf            % clear the plot window; don't use this unless you have plots -,--,:,-.
clear all      % clear all variables
format compact % condense command window output
format long g  % +, bank, hex, long, rat, short, short g, short eng

flclear fem

% COMSOL version
clear vrsn
vrsn.name = 'COMSOL 3.5';
vrsn.ext = 'a';
vrsn.major = 0;
vrsn.build = 603;
vrsn.rcs = '$Name:  $';
vrsn.date = '$Date: 2008/12/03 17:02:19 $';
fem.version = vrsn;

flbinaryfile='MP_EM_Potential_2_1_0_track_electrons.mphm';

% Geometry
clear draw
g4=flbinary('g4','draw',flbinaryfile);
draw.s.objs = {g4};
draw.s.name = {'CO1'};
draw.s.tags = {'g4'};
fem.draw = draw;
fem.geom = geomcsg(fem);
fem.mesh = flbinary('m1','mesh',flbinaryfile);

% Application mode 1
clear appl
appl.mode.class = 'Electrostatics';
appl.mode.type = 'cartesian';
appl.dim = {'V'};
appl.sdim = {'x','y','z'};
appl.name = 'es';
appl.shape = {'shlag(2,''V'')'};
appl.gporder = 4;
appl.cporder = 2;
appl.sshape = 2;
appl.border = 'off';
appl.assignsuffix = '_es';
clear prop
prop.elemdefault='Lag2';
prop.input='V';
prop.frame='ref';
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm1'};
prop.weakconstr = weakconstr;
prop.constrtype='ideal';
appl.prop = prop;
clear pnt
pnt.Q0 = 0;
pnt.V0 = 0;
pnt.type = 'Q';
pnt.ind = [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1];
appl.pnt = pnt;
clear edg
edg.Ql = 0;
edg.ind = [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, ...
  1,1,1,1,1,1];
appl.edg = edg;
clear bnd
bnd.V0 = {0,-0.1,0.1};
bnd.name = {'','Neg','Pos'};
bnd.rhos = 0;
bnd.D0 = {{0;0;0}};
bnd.epsilonrbnd = 1;
bnd.dbnd = 1;
bnd.type = {'V0','V','V'};
bnd.ind = [1,1,1,1,1,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,1];
appl.bnd = bnd;
clear equ
equ.cporder = 1;
equ.elconstrel = 'epsr';
equ.gporder = 1;
equ.Dr = {{0;0;0}};
equ.P = {{0;0;0}};
equ.epsilonr = 1;
equ.name = 'Vacuum';
equ.usage = 1;
equ.init = 0;
equ.shape = 1;
equ.rho = 0;
equ.ind = [1];
appl.equ = equ;
appl.var = {'epsilon0','8.854187817e-12'};
fem.appl{1} = appl;
fem.sdim = {'x','y','z'};
fem.frame = {'ref'};

% Simplify expressions
fem.simplify = 'on';
fem.border = 1;
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
fem.expr = {};

% Global expressions
fem.globalexpr = {};

% Functions
clear fcns
fem.functions = {};

% Solution form
fem.solform = 'weak';

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;

% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem, ...
                     'geoms',[1], ...
                     'eqvars','on', ...
                     'cplbndeq','on', ...
                     'cplbndsh','off', ...
                     'linshape',[1], ...
                     'linshapetol',0.1, ...
                     'dofversion',1);

% Retrieve solution
fem.sol=flbinary('sol1','solution',flbinaryfile);

% Save current fem structure for restart purposes
fem0=fem;

% Application mode 1
clear appl
appl.mode.class = 'Electrostatics';
appl.mode.type = 'cartesian';
appl.dim = {'V'};
appl.sdim = {'x','y','z'};
appl.name = 'es';
appl.shape = {'shlag(2,''V'')'};
appl.gporder = 4;
appl.cporder = 2;
appl.sshape = 2;
appl.border = 'off';
appl.assignsuffix = '_es';
clear prop
prop.elemdefault='Lag2';
prop.input='V';
prop.frame='ref';
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm1'};
prop.weakconstr = weakconstr;
prop.constrtype='ideal';
appl.prop = prop;
clear pnt
pnt.Q0 = 0;
pnt.V0 = 0;
pnt.type = 'Q';
pnt.ind = [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1];
appl.pnt = pnt;
clear edg
edg.Ql = 0;
edg.ind = [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, ...
  1,1,1,1,1,1];
appl.edg = edg;
clear bnd
bnd.V0 = {0,-0.1,0.1};
bnd.name = {'','Neg','Pos'};
bnd.rhos = 0;
bnd.D0 = {{0;0;0}};
bnd.epsilonrbnd = 1;
bnd.dbnd = 1;
bnd.type = {'V0','V','V'};
bnd.ind = [1,1,1,1,1,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,1];
appl.bnd = bnd;
clear equ
equ.cporder = 1;
equ.elconstrel = 'epsr';
equ.gporder = 1;
equ.Dr = {{0;0;0}};
equ.P = {{0;0;0}};
equ.epsilonr = 1;
equ.name = 'Vacuum';
equ.usage = 1;
equ.init = 0;
equ.shape = 1;
equ.rho = 0;
equ.ind = [1];
appl.equ = equ;
appl.var = {'epsilon0','8.854187817e-12'};
fem.appl{1} = appl;
fem.sdim = {'x','y','z'};
fem.frame = {'ref'};

% Simplify expressions
fem.simplify = 'on';
fem.border = 1;
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
fem.expr = {};

% Global expressions
fem.globalexpr = {};

% Functions
clear fcns
fem.functions = {};

% Solution form
fem.solform = 'weak';

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;

% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem, ...
                     'geoms',[1], ...
                     'eqvars','on', ...
                     'cplbndeq','on', ...
                     'cplbndsh','off', ...
                     'linshape',[1], ...
                     'linshapetol',0.1);

% Solve problem
fem.sol=femstatic(fem, ...
                  'method','eliminate', ...
                  'nullfun','auto', ...
                  'complexfun','off', ...
                  'matherr','on', ...
                  'solfile','off', ...
                  'conjugate','off', ...
                  'symmetric','auto', ...
                  'solcomp',{'V'}, ...
                  'outcomp',{'V'}, ...
                  'rowscale','on', ...
                  'blocksize','auto', ...
                  'reacf','on', ...
                  'ntol',1.0E-6, ...
                  'maxiter',25, ...
                  'nonlin','auto', ...
                  'damping','on', ...
                  'hnlin','off', ...
                  'linsolver','gmres', ...
                  'itol',1.0E-6, ...
                  'rhob',400.0, ...
                  'maxlinit',10000, ...
                  'itrestart',50, ...
                  'prefuntype','left', ...
                  'prefun','amg', ...
                  'prepar',{'iter',2,'mgcycle','v','mglevels',6,'maxcoarsedof',5000,'amgauto',3,'presmooth','sor','presmoothpar',{'iter',2,'relax',1.0,'sorblocked','on'},'postsmooth','soru','postsmoothpar',{'iter',2,'relax',1.0,'sorblocked','on'},'csolver','umfpack','csolverpar',{'droptol',0.0,'thresh',0.1,'umfalloc',0.7}}, ...
                  'uscale','auto', ...
                  'mcase',0);

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'V','cont','internal','recover','off','unit','V'}, ...
         'triedgestyle','none', ...
         'trifacestyle','interp', ...
         'tribar','on', ...
         'trimap','Rainbow', ...
         'trimapstyle','auto', ...
         'flowdata',{'Ex_es','Ey_es','Ez_es'}, ...
         'flowcolor',[1.0,0.0,0.0], ...
         'flowlines',30, ...
         'flowtype','line', ...
         'flowtol',0.01, ...
         'flowstattol',0.01, ...
         'flowlooptol',0.01, ...
         'flowmaxtime',Inf, ...
         'flowmaxsteps',5000, ...
         'flowback','on', ...
         'flownormal','off', ...
         'partmasstype','mass', ...
         'partplotas','lines', ...
         'parttstart','auto', ...
         'parttype','line', ...
         'partdroptimes','auto', ...
         'partstart',{[10],[10],[10]}, ...
         'partlinecolor',[0.0,0.6,1.0], ...
         'partdata',{'partq*Ex_es','partq*Ey_es','partq*Ez_es'}, ...
         'partvelstart',{'5e2','5e2','5e2'}, ...
         'partmass','9.1095e-31', ...
         'const',{'partq','-1.602e-19'}, ...
         'partrtol',0.001, ...
         'partatol','auto', ...
         'parthstart','auto', ...
         'parthmax','auto', ...
         'partstatictend','auto', ...
         'partmaxsteps','auto', ...
         'partedgetol',0.001, ...
         'partvelvar',{'partu','partv','partw'}, ...
         'parttvar','partt', ...
         'partstatic','off', ...
         'partres',5, ...
         'solnum','end', ...
         'phase',0, ...
         'title','Boundary: Electric potential [V]   Streamline: Electric field   Particle Tracing: Electromagnetic force (es)', ...
         'refine','auto', ...
         'geom','on', ...
         'geomnum',1, ...
         'bdl',{[6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21]}, ...
         'sdl',{[1]}, ...
         'complexfun','on', ...
         'matherr','off', ...
         'axisvisible','on', ...
         'axisequal','on', ...
         'grid','on', ...
         'camlight','on', ...
         'scenelight','on', ...
         'campos',[-91.31421435130811,-119.00300104368527,86.60254037844383], ...
         'camtarget',[0,0,0], ...
         'camup',[0,0,1], ...
         'camva',1.4578842528155085, ...
         'camprojection','orthographic', ...
         'transparency',1.0);

hold on;

rand ('twister', 5134879);
last_particle = 4;
for i = 1:last_particle
  particle = i;
  % Particle #1: I miss tracking the particle as it enters the POS sphere, and my track goes
  % in a different direction, though similar in shape to COMSOL's.
  % Particle #2: My beginning track is nearly identical to
  % COMSOL's, but mine leaves the POS sphere and continues, whereas COMSOL's does not.
  % Particle #3: There is good agreement at the beginning...
  % Particle #4: There is good agreement at the beginning...and I think mine makes more sense...
  % Particle #5: There is good agreement at the beginning...
  % For some particles there is good agreement, and for some there is not. I don't know why...

  % I presume that these are artifacts of the slight differences in calculations twixt COMSOL and me
  % that cause significantly different behavior near the center of the POS sphere.
  part_loc_x = [ 0.25+(rand*0.5)];
  part_loc_y = [-0.25+(rand*0.5)];
  part_loc_z = [-0.25+(rand*0.5)];

  partvelstart_x = num2str(5e2, '%10.5e\n');
  partvelstart_y = num2str(5e2, '%10.5e\n');
  partvelstart_z = num2str(5e2, '%10.5e\n');

  % Comment out the if statement if you want to see all the particles.
  if i == last_particle
    postplot(fem, ...
           'flowdata',{'Ex_es','Ey_es','Ez_es'}, ...
           'flowcolor',[1.0,0.0,0.0], ...
           'flowlines',30, ...
           'flowtype','line', ...
           'flowtol',0.01, ...
           'flowstattol',0.01, ...
           'flowlooptol',0.01, ...
           'flowmaxtime',Inf, ...
           'flowmaxsteps',5000, ...
           'flowback','on', ...
           'flownormal','off', ...
           'partmasstype','mass', ...
           'partplotas','lines', ...
           'parttstart','auto', ...
           'parttype','line', ...
           'partdroptimes','auto', ...
           'partstart',{part_loc_x,part_loc_y,part_loc_z}, ...
           'partlinecolor',[0.0,0.8,0.8], ...
           'partdata',{'-partq*Ex_es','-partq*Ey_es','-partq*Ez_es'}, ...
           'partvelstart',{partvelstart_x,partvelstart_y,partvelstart_z}, ...
           'partmass','9.1095e-31', ...
           'const',{'partq','1.602e-19'}, ...
           'partrtol',0.001, ...
           'partatol','auto', ...
           'parthstart','auto', ...
           'parthmax','auto', ...
           'partstatictend','auto', ...
           'partmaxsteps','auto', ...
           'partedgetol',0.001, ...
           'partvelvar',{'partu','partv','partw'}, ...
           'parttvar','partt', ...
           'partstatic','off', ...
           'partres',5, ...
           'solnum','end', ...
           'phase',0, ...
           'refine','auto', ...
           'geom','on', ...
           'geomnum',1, ...
           'bdl',{[3,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22]}, ...
           'sdl',{[1]}, ...
           'transparency',1.0);

  end;
end;
hold off;

% Calculations for tracking particles manaully
% Assume all calcs in 3D
% Assume all values in SI
% Note that values come from the COMSOL fem structure, and are included as hard values for
% testing...if the model is changed, these constants must be updated.
% F = q*E
% F = m*a           > a = F / m
% vf = vi + a*t
% pf = pi + vi*t + a*t*t/2 OR pf = pi + (vf + vi)*t/2

% Ex = postinterp (fem, 'Ex_es', [part_loc_x part_loc_y part_loc_z]'); % Note that the matrix is a 3D col vector
% Ey = postinterp (fem, 'Ey_es', [part_loc_x part_loc_y part_loc_z]');
% Ez = postinterp (fem, 'Ez_es', [part_loc_x part_loc_y part_loc_z]');

q  = 1.60217646e-19;    % charge of an electron
me = 9.10938188e-31;    % mass   of an electron
tstep = 1.85e-5;         % time interval twixt steps

pxyz  = [ part_loc_x part_loc_y part_loc_z ];
part_loc (1,:) = pxyz;

partvelstart_xn = str2double (partvelstart_x);
partvelstart_yn = str2double (partvelstart_y);
partvelstart_zn = str2double (partvelstart_z);
vxyzi = [ partvelstart_xn partvelstart_yn partvelstart_zn ];

for i = 2:30 % these are particles locations and begin @ 2 because the initial (above) is 1
  Ex = postinterp (fem, 'Ex_es', pxyz');
  Ey = postinterp (fem, 'Ey_es', pxyz');
  Ez = postinterp (fem, 'Ez_es', pxyz');
  Exyz = [ Ex Ey Ez ];

  % The E-Field inside the spheres may be 0...perhaps I can make them smaller to test this...
  if isnan (sum (Exyz))
    Exyz = 0.0;
  end;

  if ~isnan (sum (Exyz))
    Fxyz = -q * Exyz;
    axyz = Fxyz / me;
    vxyz = vxyzi + axyz * tstep;
    fnan = 1;
  else
    fnan = 0.1; % just a technique to slow things down when COMSOL returns a NaN E-field
  end

  pxyz = pxyz + 0.5 * (vxyz + vxyzi) * tstep * fnan;
  if particle == 0 % set to zero to disable this code block; set to 1 to agree with COMSOL #1
    % These locations are fixed to agree with COMSOL's beginning particle locations just to see if
    % the remainder of the path is similar. It seems that the E-field reported is not the same as
    % that used by COMSOL by some small difference OR there are other factors, such as adaptive
    % indexing, that are affecting calculations at each step.
    % You can see that if I follow COMSOL's track to begin with then the remainder of the track is 
    % very different. Comment out these lines, and the track makes more sense, even if it doesn't
    % agree with COMSOL.
    if i == 2
      pxyz = [ 0.593 -0.06129 0.008342 ];
    end;
    if i == 3
      pxyz = [ 0.5198 -0.01152 0.001629 ];
    end;
    if i == 4
      pxyz = [ 0.4219 0.05277 -0.007623 ];
    end;
    if i == 5
      pxyz = [ 0.3292 0.1165 -0.01258 ];
    end;
  end;
  vxyzi = vxyz;

  part_loc (i,:) = pxyz;

end

hold on;
x = part_loc (:,1);
y = part_loc (:,2);
z = part_loc (:,3);
plot3 (x,y,z,'--b');
hold off;

% Command line steps for testing:
% In this simulation,
% x is on a line twixt the charges; positive to the R
% y is the horizontal bisector of the x-line; 
% z is the vertical bisector of the x-line

% clc;
% % Starting location of the first particle - from random number generator w fixed seed
% part_loc_x =  0.619331657992472;
% part_loc_y = -0.0843199267846307;
% part_loc_z =  0.00764493003672306;
% pxyz = [ part_loc_x part_loc_y part_loc_z ]
% part_loc = pxyz;
% 
% % Starting velocity of the first particle - from random number generator w fixed seed
% partvelstart_xn = 3.12844e+004;
% partvelstart_yn = 7.30839e+004;
% partvelstart_zn = 6.44219e+004;
% vxyzi = [ partvelstart_xn partvelstart_yn partvelstart_zn ]
% 
% q  = 1.60217646e-19;    % charge of an electron
% me = 9.10938188e-31;    % mass   of an electron
% %tstep = 7.4e-8;         % time interval twixt steps
% tstep = 1.3e-7;         % time interval twixt steps
% 
% % Starting E-field @ px:py:pz - from COMSOL
% Ex =  19.9067991196741;
% Ey = -16.0452915489098;
% Ez =   0.6348407668918;
% Exyz = [ Ex Ey Ez ]
% 
% % Fx = q * Ex;
% % Fy = q * Ey;
% % Fz = q * Ez;
% Fxyz = -q * Exyz
% 
% % ax = Fx / me;
% % ay = Fy / me;
% % az = Fz / me;
% axyz = Fxyz / me
% 
% % vx = vx + ax * tstep;
% % vy = vy + ay * tstep;
% % vz = vz + az * tstep;
% vxyz = vxyzi + axyz * tstep
% 
% % px = px + vx*tstep;
% % py = py + vy*tstep;
% % pz = pz + vz*tstep;
% % pxyz = pxyz + vxyz * tstep + 0.5 * axyz * (tstep * tstep)
% pxyz = pxyz + 0.5* (vxyz + vxyzi) * tstep
% vxyzi = vxyz; % for coompleteness, because we want to do this in a loop
% 
% part_loc (2,:) = pxyz
% target_loc = [ 0.593  -0.06129  0.008342 ]
% target_loc - part_loc (2,:)
