% 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')

SimplifyMe=true;  % assumed that nspin and conductivities does not depend on the spin polariztion

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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');
%%%%%%%%%%%%%%%%%%% constant for running-wave current
model.func.create('int1', 'Interpolation');
model.func('int1').label('ADiffusion');
model.func('int1').set('source', 'file');
model.func('int1').set('filename', 'C:\myData\ADiffuse.dat');
model.func('int1').importData;
model.func('int1').set('funcname', 'ADiffusion');

model.func.create('int2', 'Interpolation');
model.func('int2').label('BInjection');
model.func('int2').set('source', 'file');
model.func('int2').set('filename', 'C:\myData\BIjection.dat');
model.func('int2').importData;
model.func('int2').set('funcname', 'BInjection');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ferro metal
sp0F=0.7;
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.8*ConductivityF');  
model.param.set('SpinConductivityF', '1*ConductivityF');  % use 1* for running-wave electron current

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

model.param.set('SpinLifeF', '30[ps]', 'spin life time');
model.param.set('nSpinXF', 'nSpin0F*(1.1428+1*0.26321)', 'the number of spin states for aprox spin diffusion length');

model.param.set('SpinLengthF', 'sqrt(SpinConductivityF*Vt*SpinLifeF/nSpinXF/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.8*ConductivityN');  
model.param.set('SpinConductivityN', '1*ConductivityN');  % % use 1* for running-wave electron current

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

model.param.set('SpinLifeN', '30[ps]', 'spin life time');
model.param.set('nSpinXN', 'nSpin0N*(1.1428+1*0.26321)', 'the number of spin states for aprox spin diffusion length');
model.param.set('SpinLengthN', 'sqrt(SpinConductivityN*Vt*SpinLifeN/nSpinXN/q)*sqrt(1-sp0N)', 'spin diffusion length');
model.param.set('Jinject0N', 'Vt*SpinConductivityN*ConductivityN/InjectionConductivityN*2/SpinLengthN', 'injection threshold current');

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

model.param.set('sp_proximity', 'sp0F*SpinLengthN/(SpinLengthN+SpinLengthF)', 'spin diffusion lengthpolarization at interface without 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', '100[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', '0[mA]','input charge current');
model.param.set('jch', 'Ich/WireArea','input scharge current density');

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

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

% spin polarization at the iterface
model.param.set('sp1', 'sp0F*SpinConductivityF/SpinLenRealF/(SpinConductivityF/SpinLenRealF+SpinConductivityN/SpinLenRealN+jch/Vt*(InjectionConductivityF/ConductivityF-InjectionConductivityN/ConductivityN))', 'spin polarization at the interface');


%%%%%%%%%%%%%%%%%%%%%%%
%%%% 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*BInjection(sp)', 'injection conductivity');
model.variable('var1').set('SpinConductivity', 'SpinConductivityF*ADiffusion(sp)', 'spin conductivity');
 if SimplifyMe 
     model.variable('var1').set('InjectionConductivity', 'InjectionConductivityF*1', 'injection conductivity');
     model.variable('var1').set('SpinConductivity', 'SpinConductivityF*1', 'spin conductivity');
 end

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

model.variable('var1').set('spA', 'sp0F+(sp1-sp0F)*exp(-x/SpinLenRealF)', 'spin polarization analitical');



% 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*BInjection(sp)', 'injection conductivity');
model.variable('var2').set('SpinConductivity', 'SpinConductivityN*ADiffusion(sp)', 'spin conductivity');
 if SimplifyMe 
     model.variable('var2').set('InjectionConductivity', 'InjectionConductivityN*1', 'injection conductivity');
     model.variable('var2').set('SpinConductivity', 'SpinConductivityN*1', 'spin conductivity');
 end

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

model.variable('var2').set('spA', 'sp0N+(sp1-sp0N)*exp(x/SpinLenRealN)', 'spin polarization analitical');

% 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');
if SimplifyMe 
    model.variable('var0').set('nSpin', 'nSpin0*(1.1428+1*0.26321)', 'it depends on spin polarization');
end

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

model.physics('c1').feature('init1').set('spV', 'spA*Vt');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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').setIndex('useDirichletCondition', '1', 1);  %spin
model.physics('c1').feature('dir1').setIndex('useDirichletCondition', '1', 0);
model.physics('c1').feature('dir1').set('r', 1, '0');   %  10 uV applied
model.physics('c1').feature('dir1').set('r', 2, 'sp0*Vt');  
model.physics('c1').feature('dir1').setIndex('useDirichletCondition', '0', 1);  % switch off
%model.physics('c1').feature('dir1').active(false);
end



if 1==1  % non-mag site
model.physics('c1').feature.create('dir2', 'DirichletBoundary', 0);
model.physics('c1').feature('dir2').selection.set([1]);
model.physics('c1').feature('dir2').setIndex('useDirichletCondition', '1', 1);  % spin
model.physics('c1').feature('dir2').setIndex('useDirichletCondition', '0', 0);
model.physics('c1').feature('dir2').set('r', 2, 'sp0*Vt');   
model.physics('c1').feature('dir2').setIndex('useDirichletCondition', '0', 1);  % switch off
%model.physics('c1').feature('dir1').active(false);
end


% ferro
if 1==1
model.physics('c1').feature.create('flux1', 'FluxBoundary', 0);
model.physics('c1').feature('flux1').selection.set([3]);
model.physics('c1').feature('flux1').set('g', 2, 'sp0*Vt/SpinLengthF*SpinConductivity*UnitC');   % transperent boundary condition 
%++sp0*jch*UnitC/Conductivity*InjectionConductivity
model.physics('c1').feature('flux1').setIndex('q', '1/SpinLengthF*SpinConductivity*UnitC+jch*UnitC/Conductivity/Vt*InjectionConductivity', 3);  % transperent boundary condition 
model.physics('c1').feature('flux1').set('g', 1, 'jch*UnitC');   %  charge flux at boundary
end



% non mag
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, 'sp0*Vt/SpinLengthN*SpinConductivity*UnitC');   % transperent boundary condition 
model.physics('c1').feature('flux2').setIndex('q', '1/SpinLengthN*SpinConductivity*UnitC-jch*UnitC/Conductivity/Vt*InjectionConductivity', 3);  % transperent boundary condition 
model.physics('c1').feature('flux2').set('g', 1, 'jch*UnitC');   %  charge flux at boundary

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





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

model.study('std1').run;
%%%%%%%%%%%%%%%%%%%%%%%%









%%%%%%%%%%%%%%%%%%  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 ' ])

SpinLenRealF=mphglobal(model,'SpinLenRealF','unit','um');
SpinLenRealN=mphglobal(model,'SpinLenRealN','unit','um');
sp1=mphglobal(model,'sp1');

disp(['analitical spin diffusion length ferro= ' num2str(SpinLenRealF) ' um non-mag= ' num2str(SpinLenRealN)  'um Spin polarization at interface= ' num2str(sp1)])
%%%%%%%%%%%%%%%%%%  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

  spInterface=sp.d1(j0);  % spin polarization at interface

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

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=1/a1;  % um


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


disp(['fitted spin diffusion length ferro= ' num2str(SpinLenF) ' um non-mag= ' num2str(SpinLenN)  'um. Spin polarization at interface= ' num2str(spInterface)])

%}

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

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


%%%%%%%%%%



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).feature.create('lngr4', 'LineGraph');
model.result(PName).feature('lngr4').set('expr', 'spA');
model.result(PName).feature('lngr4').set('xdata', 'expr');
model.result(PName).feature('lngr4').set('xdataexpr', 'x');
model.result(PName).feature('lngr4').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)=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).feature.create('lngr4', 'LineGraph');
model.result(PName).feature('lngr4').set('expr', 'log(abs(spA-sp0))');
model.result(PName).feature('lngr4').set('xdata', 'expr');
model.result(PName).feature('lngr4').set('xdataexpr', 'x');
model.result(PName).feature('lngr4').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)=-9;MyYmax(jPlot)=1;
model.result(PName).set('ymin', num2str(MyYmin(jPlot)));
model.result(PName).set('ymax', num2str(MyYmax(jPlot)));
model.result(PName).run;

%%%%%%%%%%%%%%%%%%%%%%%%%  accumulated charge
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;

%%%% conductivity
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',[ ' charge conductivity: blue; spin-diffusion conductivity:red' ' injection conductivity: green;  detection conductivity: light blue ']);
%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;

model.result(PName).feature.create('lngr2', 'LineGraph');
model.result(PName).feature('lngr2').set('expr', 'InjectionConductivity');
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]);

model.result(PName).feature.create('lngr3', 'LineGraph');
model.result(PName).feature('lngr3').set('expr', 'SpinConductivity');
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', 'DetectionConductivity');
model.result(PName).feature('lngr4').set('xdata', 'expr');
model.result(PName).feature('lngr4').set('xdataexpr', 'x');
model.result(PName).feature('lngr4').selection.set([1 2]);

%%%%% electrical field
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;

%{

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

% spin pola

    
    
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







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




%model
disp('All OK')