% that is 1D solution of spin and charge transport in a semiconductor
% model of TIA TIS assemblies
%both spin and charge currents at input
% 1 metal line
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')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
WireLength=1.5;
MyXmin=0;MyXmax=WireLength;

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

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

model.param.set('DetectionConductivity', '0.9*Conductivity'); 
model.param.set('InjectionConductivity', '0*Conductivity');  
model.param.set('SpinConductivity', '1.15*Conductivity');  % 

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

model.param.set('SpinLife', '30[ps]', 'spin life time');
model.param.set('SpinLength', 'sqrt(SpinConductivity*Vt*SpinLife/nSpin0/q)', 'spin diffusion length');


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

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


model.param.set('WireLength', [num2str(WireLength) '[um]'],'length of ferro wire');


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

%%%%%%%%%%%%%%%%%%%%%%%
%%%% 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', 'WireLength');



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




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




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



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  
model.physics('c1').feature.create('dir1', 'DirichletBoundary', 0);
model.physics('c1').feature('dir1').selection.set([1]);

model.physics('c1').feature('dir1').set('useDirichletCondition', 2, '1');
model.physics('c1').feature('dir1').set('useDirichletCondition', 1, '0');
model.physics('c1').feature('dir1').set('r', 2, 'spInput*Vt');   %  10 uV applied
%model.physics('c1').feature('dir1').active(false);
end

model.physics('c1').feature.create('dir2', 'DirichletBoundary', 0);
model.physics('c1').feature('dir2').selection.set([2]);
model.physics('c1').feature('dir2').set('useDirichletCondition', 2, '0');
model.physics('c1').feature('dir2').set('useDirichletCondition', 1, '1');
model.physics('c1').feature('dir2').set('r', 1, '0');   %  10 uV applied




if 1==1
model.physics('c1').feature.create('flux2', 'FluxBoundary', 0);
model.physics('c1').feature('flux2').selection.set([2]);
%model.physics('c1').feature('flux2').set('g', 2, 'js');   % input spin current
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

sp=mpheval(model,'sp');

M(1,:)=sp.p(:)*1000;  % nm distance
M(2,:)=sp.d1(:);

myV=mpheval(model,'uQ');
M(3,:)=myV.d1(:)*1000;  %mV

myQ=mpheval(model,'-Q/q');
M(4,:)=myQ.d1(:)*1E-6;  %cm-3

M=M';


   %{  
%}

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

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

lambda= mphglobal(model,'SpinLength','unit','nm');
disp(['spin diffusion length: ' num2str(lambda) ' nm'])
%%%%%%%%%%



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]);
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]);
MyTitle{jPlot}.fileNameSub='SpinPolarization';
model.result(PName).name(MyTitle{jPlot}.fileNameSub);


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]);
MyTitle{jPlot}.fileNameSub='LogSpinPolarization';
model.result(PName).name(MyTitle{jPlot}.fileNameSub);

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



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

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

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


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