% that is 1D solution of spin and charge transport in a semiconductor
% model of TIA TIS assemblies
%both spin and charge currents at input

clear all

import com.comsol.model.*
import com.comsol.model.util.*
model = ModelUtil.create(['Model' num2str(randi(200,1,1))]); 
model.name('metal spin transport 1D, proximity effect, assemblies');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('started')
ScanningParameter='Ich';   % Is Ich
ScanningUnit='[mA]';
MyScanningTitle=' Input charge current: ';

%MyScan=[-50:5:50];
MyScan=[50];
SaveToFile=false;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if SaveToFile
[MyFile, myPath]=uiputfile('*.*','save file','E:\delete me');
if MyFile==0;return;end
end



ScanMax=max(abs(MyScan));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
FeroWireLength=3;NonWireLength=5;
MyXmin=-NonWireLength+3;MyXmax=FeroWireLength-1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
model.param.set('q', '1.602176487e-19[C]', 'elementary charge');
model.param.set('T0', '300[K]', 'Lattice temperature');
model.param.set('k', '1.38e-23[J/K]', 'Boltzmann''s constant');
model.param.set('Vt', 'k*T0/q', 'Thermal voltage');
%%% for n0
model.param.set('epsilon0', '8.854e-12[F/m]', 'permitivity of vacuum');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ferro metal
sp0F=0.6;
model.param.set('sp0F', num2str(sp0F), 'equlibribrium spin polarization of the ferro metal');
model.param.set('DensityStateF', '2E22[1/cm^3/eV]', 'density of states of ferro');
model.param.set('nSpin0F', 'DensityStateF*k*T0', 'number of spin states without any spin accumulation');
model.param.set('ConductivityF', '2e7[S/m]', 'total conductivity of metal. Experimental fact.');  %Fe=1.12e7[S/m]  Cu=5.998e7[S/m]

model.param.set('DetectionConductivityF', '0*ConductivityF'); 
model.param.set('InjectionConductivityF', '0.9*ConductivityF');  
model.param.set('SpinConductivityF', '1.15*ConductivityF');  % 

model.param.set('epsilonF', '11.8', 'permeativity of metal ?');

model.param.set('SpinLifeF', '30[ps]', 'spin life time');
model.param.set('SpinLengthF', 'sqrt(SpinConductivityF*Vt*SpinLifeF/nSpin0F/q)*sqrt(1-sp0F)', 'spin diffusion length');
model.param.set('Jinject0F', 'Vt*SpinConductivityF*ConductivityF/InjectionConductivityF*2/SpinLengthF', 'injection threshold current');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% non-magnetic metal
sp0N=0;
model.param.set('sp0N', num2str(sp0N), 'equlibribrium spin polarization of the ferro metal');
model.param.set('DensityStateN', '2E22[1/cm^3/eV]', 'density of states of ferro');
model.param.set('nSpin0N', 'DensityStateN*k*T0', 'number of spin states without any spin accumulation');
model.param.set('ConductivityN', '2e7[S/m]', 'total conductivity of metal. Experimental fact.');  %Fe=1.12e7[S/m]  Cu=5.998e7[S/m]

model.param.set('DetectionConductivityN', '0*ConductivityN'); 
model.param.set('InjectionConductivityN', '0.9*ConductivityN');  
model.param.set('SpinConductivityN', '1.15*ConductivityN');  % 

model.param.set('epsilonN', '11.8', 'permeativity of metal ?');

model.param.set('SpinLifeN', '30[ps]', 'spin life time');

model.param.set('SpinLengthN', 'sqrt(SpinConductivityN*Vt*SpinLifeN/nSpin0N/q)*sqrt(1-sp0N)', 'spin diffusion length');
model.param.set('Jinject0N', 'Vt*SpinConductivityN*ConductivityN/InjectionConductivityN*2/SpinLengthN', 'injection threshold current');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
spInput=sp0F;
model.param.set('spInput', num2str(spInput),'input spin Polarization');

model.param.set('WireArea', '0.01[um^2]','wire cross sectional area' );


model.param.set('FeroWireLength', [num2str(FeroWireLength) '[um]'],'length of ferro wire');
model.param.set('NonWireLength', [num2str(NonWireLength) '[um]'],'length of non-magnetic metal wire');

model.param.set('WireLength', 'FeroWireLength+NonWireLength','total wire length');

model.param.set('Is', '2[uA]','input spin current');
model.param.set('js', 'Is/WireArea','input spin current density');
model.param.set('Va', '100[mV]','applied voltage');




%model.param.set('Ich', '-10[uA]','input charge current');
model.param.set('jch', 'Ich/WireArea','input scharge current density');

model.param.set('SpinLenRealN', 'SpinLengthN/(jch/Jinject0N+sqrt(1+(jch/Jinject0N)^2))', 'spin diffusion length real in non-mag');
model.param.set('SpinLenRealF', 'SpinLengthF/(-jch/Jinject0F+sqrt(1+(jch/Jinject0F)^2))', 'spin diffusion length real in ferro');

model.param.set('sp1', 'sp0F*SpinLenRealN/(SpinLenRealN+SpinLenRealF)', 'spin polarization at the interface');

model.param.set('UnitC', '1[m*s^3*A/kg]/1[S/m]', 'unit conversion');

%%%%%%%%%%%%%%%%%%%%%%%
%%%% geometry



geom1 = model.geom.create('geom1', 1);
model.geom('geom1').lengthUnit([native2unicode(hex2dec('00b5'), 'Cp1252') 'm']);  % set length unit in um

model.geom('geom1').feature.create('i1', 'Interval');
model.geom('geom1').feature('i1').set('p1', '0');
model.geom('geom1').feature('i1').set('p2', 'FeroWireLength');

model.geom('geom1').selection.create('selFerro', 'CumulativeSelection');
model.geom('geom1').feature('i1').set('contributeto', 'selFerro');

model.geom('geom1').feature.create('i2', 'Interval');
model.geom('geom1').feature('i2').set('p1', '-NonWireLength');
model.geom('geom1').feature('i2').set('p2', '0');

geom1.run;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   mesh
model.mesh.create('mesh1', 'geom1');
model.mesh('mesh1').feature.create('edg1', 'Edge');
%model.mesh('mesh1').feature('size').set('hauto', '1');
model.mesh('mesh1').feature('size').set('custom', 'on');
model.mesh('mesh1').feature('size').set('hmax', 'WireLength/2000');




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  variables

%  ferro

model.variable.create('var1');model.variable('var1').model('mod1');
model.variable('var1').selection.geom('geom1', 1);
model.variable('var1').selection.set([2]);
%model.variable('var1').selection.named('geom1_selFerro_dom');
model.variable('var1').name('VariableFerro');

model.variable('var1').set('sp0', 'sp0F', 'spin polarization');
model.variable('var1').set('nSpin0', 'nSpin0F', 'number of spin states');
model.variable('var1').set('Conductivity', 'ConductivityF', 'conductivity');
model.variable('var1').set('DetectionConductivity', 'DetectionConductivityF', 'detection conductivity');
model.variable('var1').set('InjectionConductivity', 'InjectionConductivityF', 'injection conductivity');
model.variable('var1').set('SpinConductivity', 'SpinConductivityF', 'spin conductivity');

model.variable('var1').set('epsilon', 'epsilonF', 'permeativity');
model.variable('var1').set('SpinLife', 'SpinLifeF', 'spin life time');



% non-magnetic metal
model.variable.create('var2');model.variable('var2').model('mod1');
model.variable('var2').selection.geom('geom1', 1);
model.variable('var2').name('VariableNonMag');

model.variable('var2').selection.set([1]);

model.variable('var2').set('sp0', 'sp0N', 'spin polarization');
model.variable('var2').set('nSpin0', 'nSpin0N', 'number of spin states');
model.variable('var2').set('Conductivity', 'ConductivityN', 'conductivity');

model.variable('var2').set('DetectionConductivity', 'DetectionConductivityN', 'detection conductivity');
model.variable('var2').set('InjectionConductivity', 'InjectionConductivityN', 'injection conductivity');
model.variable('var2').set('SpinConductivity', 'SpinConductivityN', 'spin conductivity');

model.variable('var2').set('epsilon', 'epsilonN', 'permeativity');
model.variable('var2').set('SpinLife', 'SpinLifeN', 'spin life time');


% both
model.variable.create('var0');model.variable('var0').model('mod1');


model.variable('var0').set('sp', 'spV/Vt', 'spin polarization');
model.variable('var0').set('nSpin', 'nSpin0*(1.1428+sp*0.26321)', 'it depends on spin polarization');


model.variable('var0').set('E', 'd(uQ,x)', 'Electrical field');
model.variable('var0').set('Q', 'd(E,x)*epsilon*epsilon0', 'accumulated charge');

model.variable('var0').set('GradSpinV', 'd(spV,x)', 'Spin Electrical field');


model.variable('var0').set('IChargeDrift', 'Conductivity*E', 'Charge drift current');
model.variable('var0').set('IChargeDifussion', 'sp*DetectionConductivity*GradSpinV', 'Charge diffusion current');
model.variable('var0').set('ICharge', 'IChargeDrift+IChargeDifussion', 'Total charge drift current');


model.variable('var0').set('ISpinDrift', 'sp*InjectionConductivity*E', 'Spin current component of drift current');
model.variable('var0').set('ISpinDifussion', 'SpinConductivity*GradSpinV', 'Diffusive spin current');
model.variable('var0').set('ISpin', 'ISpinDrift+ISpinDifussion', 'Spin current');



%%%%%%%%%%%%%%  charge-spin equation 
model.physics.create('c1', 'CoefficientFormPDE', 'geom1', {'u'});
model.physics('c1').field('dimensionless').component({'uQ' 'spV'});
model.physics('c1').prop('Units').set('DependentVariableQuantity', 1, 'electricpotential');
model.physics('c1').prop('Units').set('SourceTermQuantity', 1, 'displacement');
model.physics('c1').selection.set([1 2]);



model.physics('c1').feature('cfeq1').set('c', 2, 'InjectionConductivity*sp*UnitC' );
model.physics('c1').feature('cfeq1').set('c', 4, 'SpinConductivity*UnitC');
model.physics('c1').feature('cfeq1').set('c', 1, 'Conductivity*UnitC');
model.physics('c1').feature('cfeq1').set('c', 3, 'DetectionConductivity*sp*UnitC');


model.physics('c1').feature('cfeq1').set('f', 1, '0');
model.physics('c1').feature('cfeq1').set('f', 2, 'q*nSpin/SpinLife*sp0/(1-sp0)*Vt*UnitC/Vt');  % clasical=0
model.physics('c1').feature('cfeq1').set('a', 4, 'q*nSpin/SpinLife*(1+sp0/(1-sp0))*UnitC/Vt');  % clasical Conduct/SpinLen^2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if 1==1  % ferro site
model.physics('c1').feature.create('dir1', 'DirichletBoundary', 0);
model.physics('c1').feature('dir1').selection.set([3]);
model.physics('c1').feature('dir1').set('useDirichletCondition', 2, '0');
model.physics('c1').feature('dir1').set('useDirichletCondition', 2, '1');
model.physics('c1').feature('dir1').set('r', 1, '0');   %  10 uV applied
model.physics('c1').feature('dir1').set('r', 2, 'spInput*Vt');   %  10 uV applied
%model.physics('c1').feature('dir1').active(false);
end

if 1==0
model.physics('c1').feature.create('dir2', 'DirichletBoundary', 0);
model.physics('c1').feature('dir2').selection.set([1]);
model.physics('c1').feature('dir2').set('useDirichletCondition', 2, '1[mV]');  % does not apply to spin
model.physics('c1').feature('dir2').set('r', 1, '0');   %  for charge .ground
end

if 1==0
model.physics('c1').feature.create('flux1', 'FluxBoundary', 0);
model.physics('c1').feature('flux1').selection.set([2]);
model.physics('c1').feature('flux1').set('g', 2, '-uS/SpinLenEff*Conduct');   % open boundary
model.physics('c1').feature('flux1').set('g', 1, '0');   % zero charge flux at boundary
end

if 1==1
model.physics('c1').feature.create('flux2', 'FluxBoundary', 0);
model.physics('c1').feature('flux2').selection.set([1]);
%model.physics('c1').feature('flux2').set('g', 2, 'js');   % input spin current
model.physics('c1').feature('flux2').set('g', 1, 'jch');   %  charge flux at boundary
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%





%%%%%%%%%%%%%%%%%%%%%%%
model.study.create('std1');
model.study('std1').feature.create('stat', 'Stationary');


%%%%%%%%%%%%%%%%%%%%%%%%

for j11=1:length(MyScan)
    %j1=j11;  % positive
    j1=length(MyScan)+1-j11;  % negative
model.param.set(ScanningParameter, [num2str(MyScan(j11)) ScanningUnit]);

model.study('std1').run;

InputCurrent(j1)=MyScan(j11);





%%%%%%%%%%%%%%%%%%  read spin polarization

IthF=mphglobal(model,'Jinject0F','unit','mA/um^2');
IthN=mphglobal(model,'Jinject0N','unit','mA/um^2');
disp(['threshold current ferro= ' num2str(IthF) ' mA/um^2 non-mag= ' num2str(IthN)  'mA/um^2'])
my_jch=mphglobal(model,'jch','unit','mA/um^2');
disp(['current used=  ' num2str(my_jch) ' mA/um^2 ' ])

lambdaRealF=mphglobal(model,'SpinLenRealF','unit','um');
lambdaRealN=mphglobal(model,'SpinLenRealN','unit','um');
disp(['real spin diffusion length ferro= ' num2str(lambdaRealF) ' um non-mag= ' num2str(lambdaRealN)  'um'])

%%%%%%%%%%%%%%%%%%  read spin polarization

sp=mpheval(model,'sp');
MyL=length(sp.p);
%
for j=1:length(sp.p)
    if sp.p(j)>=0
        j0=j;break;
    end
end

   

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

jNm=fix(j0*0.1);
toFitXF=sp.p(5:j0);
toFitYF=log(abs(sp.d1(jNm:j0)-sp0N))-log(abs(sp.d1(j0)-sp0N));


toFitXN=sp.p(j0:(MyL-5));
toFitYN=log(abs(sp.d1(j0:(MyL-5))-sp0F))-log(abs(sp.d1(j0)-sp0F));
%hh=figure(jPlot);plot(toFitXF,toFitYF,toFitXN,toFitYN);grid on;
%F(j1,1)=getframe(hh);
[a1 , ~, ~]=LinearFit(toFitYF,toFitXF);
SpinLenN(j1)=1/a1;  % um
SpinLen(j1,1)=MyScan(j11);
SpinLen(j1,2)=SpinLenN(j1);


[a1,~,~]=LinearFit(toFitYN,toFitXN);
SpinLenF(j1)=-1/a1;  % um
SpinLen(j1,3)=SpinLenF(j1);

disp(['fitted spin diffusion length ferro= ' num2str(SpinLenF(j1)) ' um non-mag= ' num2str(SpinLenN(j1))  'um'])
   %{  
%}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  plots

jPlot=0;
%MyChangeTitle=[ScanningParameter ' ' num2str(MyScan(j11))  ' ' ScanningUnit   ];
MyChangeTitle=[MyScanningTitle ' ' num2str(MyScan(j11))  ' ' ScanningUnit   ];


%%%%%%%%%%

if j11==1

jPlot=jPlot+1;PName=['pg' num2str(jPlot)];MyName{jPlot}=PName;
model.result.create(PName, 'PlotGroup1D');
model.result(PName).feature.create('lngr1', 'LineGraph');
model.result(PName).feature('lngr1').set('expr', 'uQ');
model.result(PName).feature('lngr1').set('xdata', 'expr');
model.result(PName).feature('lngr1').set('xdataexpr', 'x');
model.result(PName).feature('lngr1').set('unit', 'mV');
MyTitle{jPlot}.myName='Voltage  mV';
model.result(PName).set('title',[ MyChangeTitle]);
model.result(PName).set('ylabelactive', 'on');
model.result(PName).set('ylabel', MyTitle{jPlot}.myName);
model.result(PName).feature('lngr1').selection.set([1 2]);
MyTitle{jPlot}.fileNameSub='Voltage';
model.result(PName).name(MyTitle{jPlot}.fileNameSub);
model.result(PName).run;


jPlot=jPlot+1;PName=['pg' num2str(jPlot)];MyName{jPlot}=PName;
model.result.create(PName, 'PlotGroup1D');
model.result(PName).feature.create('lngr1', 'LineGraph');
model.result(PName).feature('lngr1').set('expr', 'sp');
model.result(PName).feature('lngr1').set('xdata', 'expr');
model.result(PName).feature('lngr1').set('xdataexpr', 'x');
MyTitle{jPlot}.myName='Spin Polarization';
model.result(PName).set('title',[ MyChangeTitle]); 
model.result(PName).set('ylabelactive', 'on');
model.result(PName).set('ylabel', MyTitle{jPlot}.myName);
model.result(PName).feature('lngr1').selection.set([1 2]);
MyTitle{jPlot}.fileNameSub='SpinPolarization';
model.result(PName).name(MyTitle{jPlot}.fileNameSub);

model.result(PName).feature.create('lngr2', 'LineGraph');
model.result(PName).feature('lngr2').set('expr', 'sp1*exp(x/SpinLenRealN)');
model.result(PName).feature('lngr2').set('xdata', 'expr');
model.result(PName).feature('lngr2').set('xdataexpr', 'x');
model.result(PName).feature('lngr2').selection.set([1]);

model.result(PName).feature.create('lngr3', 'LineGraph');
model.result(PName).feature('lngr3').set('expr', 'sp0F+(sp1-sp0F)*exp(-x/SpinLenRealF)');
model.result(PName).feature('lngr3').set('xdata', 'expr');
model.result(PName).feature('lngr3').set('xdataexpr', 'x');
model.result(PName).feature('lngr3').selection.set([2]);


model.result(PName).set('axislimits', 'on');
model.result(PName).set('xmin', num2str(MyXmin));
model.result(PName).set('xmax', num2str(MyXmax));
MyYmin(jPlot)=0;MyYmax(jPlot)=1;
model.result(PName).set('ymin', num2str(MyYmin(jPlot)));
model.result(PName).set('ymax', num2str(MyYmax(jPlot)));
model.result(PName).run;

jPlot=jPlot+1;PName=['pg' num2str(jPlot)];MyName{jPlot}=PName;
model.result.create(PName, 'PlotGroup1D');
model.result(PName).feature.create('lngr1', 'LineGraph');
model.result(PName).feature('lngr1').set('expr', 'log(abs(sp-sp0))');
model.result(PName).feature('lngr1').set('xdata', 'expr');
model.result(PName).feature('lngr1').set('xdataexpr', 'x');
MyTitle{jPlot}.myName='Log of Spin Polarization';
model.result(PName).set('title',[ MyChangeTitle]);
model.result(PName).set('ylabelactive', 'on'); 
model.result(PName).set('ylabel', MyTitle{jPlot}.myName); 
model.result(PName).feature('lngr1').selection.set([1 2]);
MyTitle{jPlot}.fileNameSub='LogSpinPolarization';
model.result(PName).name(MyTitle{jPlot}.fileNameSub);

model.result(PName).feature.create('lngr2', 'LineGraph');
model.result(PName).feature('lngr2').set('expr', 'log(sp1)+x/SpinLenRealN');
model.result(PName).feature('lngr2').set('xdata', 'expr');
model.result(PName).feature('lngr2').set('xdataexpr', 'x');
model.result(PName).feature('lngr2').selection.set([1]);

model.result(PName).feature.create('lngr3', 'LineGraph');
model.result(PName).feature('lngr3').set('expr', 'log(sp1-sp0F)-(x/SpinLenRealF)');
model.result(PName).feature('lngr3').set('xdata', 'expr');
model.result(PName).feature('lngr3').set('xdataexpr', 'x');
model.result(PName).feature('lngr3').selection.set([2]);

model.result(PName).set('axislimits', 'on');
model.result(PName).set('xmin', num2str(MyXmin));
model.result(PName).set('xmax', num2str(MyXmax));
MyYmin(jPlot)=-9;MyYmax(jPlot)=1;
model.result(PName).set('ymin', num2str(MyYmin(jPlot)));
model.result(PName).set('ymax', num2str(MyYmax(jPlot)));
model.result(PName).run;


jPlot=jPlot+1;PName=['pg' num2str(jPlot)];MyName{jPlot}=PName;
model.result.create(PName, 'PlotGroup1D');
model.result(PName).feature.create('lngr1', 'LineGraph');
model.result(PName).feature('lngr1').set('expr', '-Q/q');
model.result(PName).feature('lngr1').set('xdata', 'expr');
model.result(PName).feature('lngr1').set('xdataexpr', 'x');
model.result(PName).feature('lngr1').set('unit', '1/cm^3');
MyTitle{jPlot}.myName='Accumulated electrons';
model.result(PName).set('title',[ MyChangeTitle]);
model.result(PName).set('ylabelactive', 'on');
model.result(PName).set('ylabel', MyTitle{jPlot}.myName);
model.result(PName).feature('lngr1').selection.set([1 2]);
%model.result(PName).set('ylog', 'on');
MyTitle{jPlot}.fileNameSub='AccumulatedCharge';
model.result(PName).name(MyTitle{jPlot}.fileNameSub);
model.result(PName).run;


jPlot=jPlot+1;PName=['pg' num2str(jPlot)];MyName{jPlot}=PName;
model.result.create(PName, 'PlotGroup1D');
model.result(PName).feature.create('lngr1', 'LineGraph');
model.result(PName).feature('lngr1').set('expr', 'Conductivity');
model.result(PName).feature('lngr1').set('xdata', 'expr');
model.result(PName).feature('lngr1').set('xdataexpr', 'x');
MyTitle{jPlot}.myName='Conductivity';
model.result(PName).set('title',[ MyChangeTitle]);
model.result(PName).set('ylabelactive', 'on');
model.result(PName).set('ylabel', MyTitle{jPlot}.myName);
model.result(PName).feature('lngr1').selection.set([1 2]);
MyTitle{jPlot}.fileNameSub='Conductivity';
model.result(PName).name(MyTitle{jPlot}.fileNameSub);
model.result(PName).run;

jPlot=jPlot+1;PName=['pg' num2str(jPlot)];MyName{jPlot}=PName;
model.result.create(PName, 'PlotGroup1D');
model.result(PName).feature.create('lngr1', 'LineGraph');
model.result(PName).feature('lngr1').set('expr', 'E');
model.result(PName).feature('lngr1').set('xdata', 'expr');
model.result(PName).feature('lngr1').set('xdataexpr', 'x');
MyTitle{jPlot}.myName='Electrical field';
model.result(PName).set('title',[ MyChangeTitle]);
model.result(PName).set('ylabelactive', 'on');
model.result(PName).set('ylabel', MyTitle{jPlot}.myName);
model.result(PName).feature('lngr1').selection.set([1 2]);
MyTitle{jPlot}.fileNameSub='ElectricalField';
model.result(PName).name(MyTitle{jPlot}.fileNameSub);
model.result(PName).run;

jPlot=jPlot+1;PName=['pg' num2str(jPlot)];MyName{jPlot}=PName;
model.result.create(PName, 'PlotGroup1D');
model.result(PName).feature.create('lngr1', 'LineGraph');
model.result(PName).feature('lngr1').set('expr', 'WireArea*ICharge');
model.result(PName).feature('lngr1').set('xdata', 'expr');
model.result(PName).feature('lngr1').set('xdataexpr', 'x');
MyTitle{jPlot}.myName='Charge current in mA';
model.result(PName).feature('lngr1').set('unit','mA');
model.result(PName).set('title',[ MyChangeTitle]);
model.result(PName).set('ylabelactive', 'on');
model.result(PName).set('ylabel', MyTitle{jPlot}.myName);
model.result(PName).feature('lngr1').selection.set([1 2]);
MyTitle{jPlot}.fileNameSub='ChargeCurrent';
model.result(PName).name(MyTitle{jPlot}.fileNameSub);
model.result(PName).run;



jPlot=jPlot+1;PName=['pg' num2str(jPlot)];MyName{jPlot}=PName;
model.result.create(PName, 'PlotGroup1D');
model.result(PName).feature.create('lngr1', 'LineGraph');
model.result(PName).feature('lngr1').set('expr', '-WireArea*ISpin');
model.result(PName).feature('lngr1').set('xdata', 'expr');
model.result(PName).feature('lngr1').set('xdataexpr', 'x');
MyTitle{jPlot}.myName=['Spin current in mA total(blue); drift(green) diffuse(red)'];
%MyTitle{jPlot}.myName=['Spin current in mA ; total(blue); drift(green) ' num2str(InputISpinDrift*1E6) 'uA; diffuse(red): ' num2str(InputISpinDiffus*1E6) ' uA']);
model.result(PName).feature('lngr1').set('unit','mA');
model.result(PName).set('title',[ MyChangeTitle]);
model.result(PName).set('ylabelactive', 'on');
model.result(PName).set('ylabel', MyTitle{jPlot}.myName);
model.result(PName).feature('lngr1').selection.set([1 2]);
MyTitle{jPlot}.fileNameSub='SpinCurrent';
model.result(PName).name(MyTitle{jPlot}.fileNameSub);



model.result(PName).feature.create('lngr2', 'LineGraph');
model.result(PName).feature('lngr2').set('expr', '-WireArea*ISpinDrift');
model.result(PName).feature('lngr2').set('xdata', 'expr');
model.result(PName).feature('lngr2').set('xdataexpr', 'x');
model.result(PName).feature('lngr2').set('unit','mA');
model.result(PName).feature('lngr2').selection.set([1 2]);



model.result(PName).feature.create('lngr3', 'LineGraph');
model.result(PName).feature('lngr3').set('expr', '-WireArea*ISpinDifussion');
model.result(PName).feature('lngr3').set('xdata', 'expr');
model.result(PName).feature('lngr3').set('xdataexpr', 'x');
model.result(PName).feature('lngr3').set('unit','mA');
model.result(PName).feature('lngr3').selection.set([1 2]);

model.result(PName).set('axislimits', 'on');
model.result(PName).set('xmin', num2str(MyXmin));
model.result(PName).set('xmax', num2str(MyXmax));
%MyYmin(jPlot)=-1.2;MyYmax(jPlot)=0.8;  % current -10  10
MyYmin(jPlot)=-30;MyYmax(jPlot)=30;
model.result(PName).set('ymin', num2str(MyYmin(jPlot)));
model.result(PName).set('ymax', num2str(MyYmax(jPlot)));

%{
model.result(PName).feature.create('lngr4', 'LineGraph');
model.result(PName).feature('lngr4').set('expr', 'WireArea*js*exp(-x/SpinLen)');
model.result(PName).feature('lngr4').set('xdata', 'expr');
model.result(PName).feature('lngr4').set('xdataexpr', 'x');
model.result(PName).feature('lngr4').set('unit','uA');
model.result(PName).feature('lngr4').selection.set([1]);
%}
model.result(PName).run;
%{
jPlot=jPlot+1;PName=['pg' num2str(jPlot)];MyName{jPlot}=PName;
model.result.create(PName, 'PlotGroup1D');
model.result(PName).feature.create('lngr1', 'LineGraph');
model.result(PName).feature('lngr1').set('expr', 'log10(abs(WireArea*ISpin/1[uA]))');
model.result(PName).feature('lngr1').set('xdata', 'expr');
model.result(PName).feature('lngr1').set('xdataexpr', 'x');
MyTitle{jPlot}.myName='log10(log10(SpinCurrent/1uA); total-blue; drift-green; diffuse-red';
%model.result(PName).feature('lngr1').set('unit','mA');
model.result(PName).feature('lngr1').selection.set([1 2]);
model.result(PName).name('Log of Spin current');

% spin component of drift current
model.result(PName).feature.create('lngr2', 'LineGraph');
model.result(PName).feature('lngr2').set('expr', 'log10(WireArea*ISpinDrift/1[uA])');
model.result(PName).feature('lngr2').set('xdata', 'expr');
model.result(PName).feature('lngr2').set('xdataexpr', 'x');
model.result(PName).feature('lngr2').selection.set([1 2]);

% diffusive spin current
model.result(PName).feature.create('lngr3', 'LineGraph');
model.result(PName).feature('lngr3').set('expr', 'log10(WireArea*ISpinDifussion/1[uA])');
model.result(PName).feature('lngr3').set('xdata', 'expr');
model.result(PName).feature('lngr3').set('xdataexpr', 'x');
model.result(PName).feature('lngr3').selection.set([1 2]);


model.result(PName).feature.create('lngr4', 'LineGraph');
model.result(PName).feature('lngr4').set('expr', 'log10(WireArea*js/1[uA])-x/SpinLenEff*log10(exp(1))');
model.result(PName).feature('lngr4').set('xdata', 'expr');
model.result(PName).feature('lngr4').set('xdataexpr', 'x');
model.result(PName).feature('lngr4').selection.set([1]);

%}

%model.result(PName).run;
PlotNumber=jPlot;
end


for j=1:PlotNumber
h(j)=figure(j); mphplot(model,MyName{j});
model.result(MyName{j}).set('title',[ MyChangeTitle]);

% spin pola
if j==2
    ylabel('Spin Polarization, sp' ,'FontSize', 14)
    hold on
    rectangle('Position',[MyXmin 0.92 -MyXmin 0.08],'FaceColor',[1 1 0.8]); %rectangle('Position',[x,y,w,h])
    rectangle('Position',[0 0.92 MyXmax 0.08],'FaceColor',[1 0.9 1]);
    text(MyXmin*0.9,0.95,' non-magnetic ','FontSize',16)
    text(0,0.95,' ferro','FontSize',16)
    
    k=abs(MyScan(j11)/ScanMax);
    MyXLen=MyXmax+abs(MyXmin);
    ymax=1;
    
    if MyScan(j11)<0
 arrows(0,0.85*ymax,k*MyXLen*0.3,-90,k*[ymax*0.03,0.4,0.4,ymax*0.01],'FaceColor','g');
    else
        arrows(0,0.85*ymax,k*MyXLen*0.3,90,k*[ymax*0.03,0.4,0.4,ymax*0.01],'FaceColor','r');
    end
    
    hold off
    
    
end


if j==3 || j==8
    hold on
    ymax=MyYmax(j);ymin=MyYmin(j); dy=ymax-ymin;
    rectangle('Position',[MyXmin ymin+dy*0.92 -MyXmin dy*0.08],'FaceColor',[1 1 0.8]); %rectangle('Position',[x,y,w,h])
    rectangle('Position',[0 ymin+dy*0.92 MyXmax dy*0.08],'FaceColor',[1 0.9 1]);
    text(MyXmin*0.9,ymin+dy*0.95,' non-magnetic ','FontSize',16)
    text(0,ymin+dy*0.95,' ferro','FontSize',16)
    
    k=abs(MyScan(j11)/ScanMax);
    MyXLen=MyXmax+abs(MyXmin);
    
    
    if MyScan(j11)<0
 arrows(0,ymin+dy*0.85,k*MyXLen*0.3,-90,k*[dy*0.03,0.4,0.4,dy*0.01],'FaceColor','g');
    else
        arrows(0,ymin+dy*0.85,k*MyXLen*0.3,90,k*[dy*0.03,0.4,0.4,dy*0.01],'FaceColor','r');
    end
    
    hold off
    
end

if j==3 
 ylabel('Log of Spin Polarization','FontSize',15)  ;
end
if j==8 
 ylabel('Spin current, mA','FontSize',15)  ;
end

xlabel('coordinate across contact, \mum','FontSize',14)
title(MyChangeTitle,'FontSize',16)

if j==8 
    title({MyChangeTitle; 'total(blue); drift(green) diffuse(red)'},'FontSize',14)

end

F(j1,j)=getframe(h(j));

end


end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if SaveToFile
    %close all;
    if 1==1 % no picture
    for j2=1:j
        F1(:,1)=F(:,j2);
    save([myPath MyFile MyTitle{j2}.fileNameSub 'F1' ],'F1');  %save images frame
  clear F1;
    end
    end
end

  %{  
   
 fid = fopen([myPath MyFile '.txt' ],'w');
    fprintf(fid,['InputCurrent InputISpinDrift InputISpinDiffus DriftRatio DiffusRatio SpinLenEff SpinLenDrift SpinLenDiffus SpinLenSpinV \n']);
    
    for j2=1:length(InputCurrent)        
 fprintf(fid,'%0.5g  %0.5g    %0.5g  %0.5g  %0.5g  %0.5g  %0.5g %0.5g %0.5g \n',InputCurrent(j2), InputISpinDrift(j2), InputISpinDiffus(j2), DriftRatio(j2), DiffusRatio(j2), SpinLenEff(j2), SpinLen1(j2), SpinLen2(j2), SpinLen3(j2));
    end
    
 fclose(fid); 
 
end
%}
    
h1=figure(1);movie(h1,F(:,2),100,3);
close all;
disp('All OK')

%model
disp('All OK')


%%%%%%% get data