% that is 2D solution of spin excited in metal (Crooker experiment)
% both spin and charge currents at input

clear all


import com.comsol.model.*
import com.comsol.model.util.*
model = ModelUtil.create(['Model' num2str(randi(900,1,1))]); 
model.name('Crooker metal');
%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('started')
ScanningParameter='Va';   % Is Ich
ScanningUnit='[mV]';
MyScanningTitle=' Applied voltage: ';

MyScan=[-2000:200:2000];
SaveToFile=true;

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



ScanMax=max(abs(MyScan));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
GaussianBeam='y';


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%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');
model.param.set('epsilon0', '8.854e-12[F/m]', 'permitivity of vacuum');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% non-magnetic metal

model.param.set('DensityStateN', '2E22[1/cm^3/eV]', 'density of states of ferro');
model.param.set('nSpin0', 'DensityStateN*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*Conductivity'); 
model.param.set('InjectionConductivity', '0.9*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');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% p



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

model.param.set('CalculationLength', 'SpinLength*20'  );
model.param.set('CalculationWidth', 'SpinLength*20'  );


model.param.set('FilmThickness', '1[um]','Film thickness');


model.param.set('Va', '-2000[mV]','applied voltage');


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

model.param.set('BeamDiameter', 'SpinLength' );
model.param.set('TimePump', '40[ps]' );
model.param.set('sp0Pump', 'SpinLife/(SpinLife+TimePump)' );


% for semicond
%{
model.param.set('BeamPower', '5[uW]', 'Laser Beam Power');
model.param.set('BeamLambda', '800[nm]', 'Laser Wavelength');
model.param.set('PhotonEnergy', 'c_const*h_const/BeamLambda','Single Photon Energy');
model.param.set('PhotonTime', 'BeamPower/PhotonEnergy');
model.param.set('SpinPhotonEfficiency', '0.34','Photon to spin conversion efficiency');  
% 50 when light hole is not splitted  % 32 normal refractivity meaning 68 %
% go through   0.5*0.68=0.34%

model.param.set('GenSpinTime', 'PhotonTime*SpinPhotonEfficiency','number of spin generated per second');
%}

%%%%%%%%%%%%%%  gausian beam
if GaussianBeam=='y'

model.func.create('an1', 'Analytic');
model.func('an1').set('funcname', 'fGauss');
model.func('an1').set('args', 'x, y');
model.func('an1').set('expr', 'exp(-(x^2+y^2)/BeamDiameter^2)');  % that is normalized gaussian Laser2 is square of R

end
%%%%%%%%%%%%%%%%%%%%%%%
%%%% geometry

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

%%%%%%%%%%%% calculation window
r1=model.geom('geom1').feature.create('r1', 'Rectangle');
r1.set('size',{'CalculationLength' 'CalculationWidth'});
r1.set('pos',{'-CalculationLength/2' '-CalculationWidth/2'});
geom1.run
model.geom('geom1').feature().create('sel1', 'Selection');
model.geom('geom1').feature('sel1').selection('selection').set('r1',1);
geom1.run


model.geom('geom1').feature.create('c1', 'Circle');
model.geom('geom1').feature('c1').setIndex('pos', '0', 0);  % x of center
model.geom('geom1').feature('c1').setIndex('pos', '0', 1);   % y-of center
if GaussianBeam=='y'
 model.geom('geom1').feature('c1').set('r', '3*BeamDiameter/2');   % in case of gaussian calculation in 3 R
else
model.geom('geom1').feature('c1').set('r', 'BeamDiameter/2');
end

model.geom('geom1').runAll;
model.geom('geom1').feature().create('sel2', 'Selection');
model.geom('geom1').feature('sel2').selection('selection').set('c1',1);
geom1.run

%%%%%%%%%%%%%%%%%%
Main_Sel=model.selection('geom1_sel1').entities(2);
Beam_Sel=model.selection('geom1_sel2').entities(2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   mesh
%%%%%%%
% create mesh
%  create mesh
mesh1=model.mesh.create('mesh1', 'geom1');
meshSize=mesh1.feature('size');
meshSize.set('hauto', '6');
mesh1.feature.create('ftri1', 'FreeTri');

%model.mesh('mesh1').feature('size').set('custom', 'on');
%model.mesh('mesh1').feature('size').set('hmax', 'SpinLen/150');

mesh1.run;
mesh1.data.transferMesh;

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


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

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

model.variable('var1').set('Ex', 'd(uQ,x)', 'Electrical field');
model.variable('var1').set('Ey', 'd(uQ,y)', 'Electrical field');

model.variable('var1').set('Q', '(d(Ex,x)+d(Ey,y))*epsilon*epsilon0', 'accumulated charge');

model.variable('var1').set('GradSpinVx', 'd(spV,x)', 'Spin Electrical field');
model.variable('var1').set('GradSpinVy', 'd(spV,y)', 'Spin Electrical field');



model.variable('var1').set('IChargeDriftx', 'Conductivity*Ex', 'Charge drift current');
model.variable('var1').set('IChargeDrifty', 'Conductivity*Ey', 'Charge drift current');

model.variable('var1').set('IChargeDifussionx', 'sp*DetectionConductivity*GradSpinVx', 'Charge diffusion current');
model.variable('var1').set('IChargeDifussiony', 'sp*DetectionConductivity*GradSpinVy', 'Charge diffusion current');

model.variable('var1').set('IChargex', 'IChargeDriftx+IChargeDifussionx', 'Total charge drift current');
model.variable('var1').set('IChargey', 'IChargeDrifty+IChargeDifussiony', 'Total charge drift current');



model.variable('var1').set('ISpinDriftx', 'sp*InjectionConductivity*Ex', 'Spin current component of drift current');
model.variable('var1').set('ISpinDrifty', 'sp*InjectionConductivity*Ey', 'Spin current component of drift current');

model.variable('var1').set('ISpinDifussionx', 'SpinConductivity*GradSpinVx', 'Diffusive spin current');
model.variable('var1').set('ISpinDifussiony', 'SpinConductivity*GradSpinVy', 'Diffusive spin current');

model.variable('var1').set('ISpinx', 'ISpinDriftx+ISpinDifussionx', 'Spin current');
model.variable('var1').set('ISpiny', 'ISpinDrifty+ISpinDifussiony', 'Spin current');






%%%% main varible
model.variable.create('var2');model.variable('var2').model('mod1');
model.variable('var2').selection.geom('geom1', 2);
model.variable('var2').selection.set([1]);
model.variable('var2').set('sp0', '0');

%%%% main varible
model.variable.create('var3');model.variable('var3').model('mod1');
model.variable('var3').selection.geom('geom1', 2);

model.variable('var3').selection.set([2]);
if GaussianBeam=='y'
    model.variable('var3').set('sp0', 'sp0Pump*fGauss(x,y)');   
else
model.variable('var3').set('sp0', 'sp0Pump');
end




%%%%%%%%%%%%%%  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' '0' '0' 'InjectionConductivity*sp*UnitC'});
model.physics('c1').feature('cfeq1').set('c', 4, {'SpinConductivity*UnitC' '0' '0' 'SpinConductivity*UnitC'});
model.physics('c1').feature('cfeq1').set('c', 1, {'Conductivity*UnitC' '0' '0' 'Conductivity*UnitC'});
model.physics('c1').feature('cfeq1').set('c', 3, {'DetectionConductivity*sp*UnitC' '0' '0' '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', 1);
model.physics('c1').feature('dir1').selection.set(4);
model.physics('c1').feature('dir1').set('useDirichletCondition', 2, '0');  % does not apply to spin
model.physics('c1').feature('dir1').set('r', 1, 'Va');   %  10 uV applied
%model.physics('c1').feature('dir1').active(false);

end

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




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

MyVoltage(j1)=MyScan(j11);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%result
myE = mphinterp(model,'Ex','coord',[0;0],'unit','V/cm');
myI = mphinterp(model,'IChargeDriftx','coord',[0;0],'unit','mA/um^2');
%commonTitle=[' Electrical Field: ' num2str(myE) ' V/cm Current ' num2str(myI) ' uA/um2']; 
%commonTitle=' ';
commonTitle=[' Electrical Field: ' num2str(round(myE/10)*10) ' V/cm ']; 


j=0;
j=j+1; pg=['pg' +num2str(j)];
if j11==1
model.result.create(pg, 2);
model.result(pg).feature.create('surf1', 'Surface');
model.result(pg).feature('surf1').set('expr', 'sp');
model.result(pg).set('titleactive','on');
%model.result(pg).set('title',['spin polarization ' commonTitle]);
model.result(pg).set('title',['  ']);
model.result(pg).name('SpinPolarization');
model.result(pg).feature('surf1').set('rangecoloractive', 'on');
model.result(pg).feature('surf1').set('rangecolormin', '0');
model.result(pg).feature('surf1').set('rangecolormax', '0.18');

end



model.result(pg).run;

h(j)=figure(j); mphplot(model,pg,'rangenum',1);
axis([-2.5,2.5,-2.5,2.5])
xlabel('x-coordinate, \mum','FontSize', 18)
ylabel('y-coordinate, \mum','FontSize', 18)
title({['Spin Polarization']; commonTitle}, 'FontSize',16)
F(j11,j)=getframe(h(j));
MyTitle{j}='SpinPol';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%{
j=j+1; pg=['pg' +num2str(j)];
model.result.create(pg, 2);
%model.result(pg).set('data', 'dset1');
model.result(pg).feature.create('surf1', 'Surface');
model.result(pg).feature('surf1').set('expr', 'uQ');
model.result(pg).feature('surf1').set('unit', 'uV');
model.result(pg).set('titleactive','on');
model.result(pg).set('title',[' Voltage [uV] ']);
model.result(pg).name('Volatage');
model.result(pg).run;
%figure(j); mphplot(model,pg,'rangenum',1);



j=j+1; pg=['pg' +num2str(j)];
model.result.create(pg, 2);
%model.result(pg).set('data', 'dset1');
model.result(pg).feature.create('surf1', 'Surface');
model.result(pg).feature('surf1').set('expr', 'sqrt(Ex^2+Ey^2)');
model.result(pg).feature('surf1').set('unit', 'V/cm');
model.result(pg).set('titleactive','on');
model.result(pg).set('title',[' Electrical Field [V/cm] ']);
model.result(pg).name('Field');
model.result(pg).run;
%figure(j); mphplot(model,pg,'rangenum',1);
%}
j=j+1; pg=['pg' +num2str(j)];
if j11==1
model.result.create(pg, 2);
%model.result(pg).set('data', 'dset1');
model.result(pg).feature.create('surf1', 'Surface');
model.result(pg).feature('surf1').set('expr', 'sqrt(ISpinDifussionx^2+ISpinDifussiony^2)');
model.result(pg).feature('surf1').set('unit', 'mA/um^2');
model.result(pg).set('titleactive','on');
model.result(pg).set('title',[' Diffusion Spin Current, , A/um^2   ']);
model.result(pg).name('DiffussionCurrent');

model.result(pg).feature('surf1').set('rangecoloractive', 'on');
model.result(pg).feature('surf1').set('rangecolormin', '0');
model.result(pg).feature('surf1').set('rangecolormax', '180');

model.result(pg).feature.create('arws1', 'ArrowSurface');
model.result(pg).feature('arws1').setIndex('expr', '-ISpinDifussionx', 0);
model.result(pg).feature('arws1').setIndex('expr', '-ISpinDifussiony', 1);
model.result(pg).feature('arws1').set('arrowlength', 'logarithmic');


end



model.result(pg).run;

h(j)=figure(j); mphplot(model,pg,'rangenum',1);
axis([-2.5,2.5,-2.5,2.5])
xlabel('x-coordinate, \mum','FontSize', 18)
ylabel('y-coordinate, \mum','FontSize', 18)
title({['Diffusion Spin Current, , mA/um^2']; commonTitle}, 'FontSize',16)
F(j11,j)=getframe(h(j));
MyTitle{j}='SpinDiffus';

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

j=j+1; pg=['pg' +num2str(j)];
if j11==1
model.result.create(pg, 2);
%model.result(pg).set('data', 'dset1');
model.result(pg).feature.create('surf1', 'Surface');
model.result(pg).feature('surf1').set('expr', 'sqrt(ISpinDriftx^2+ISpinDrifty^2)');
model.result(pg).feature('surf1').set('unit', 'mA/um^2');
model.result(pg).set('titleactive','on');
model.result(pg).set('title',[' Drift Spin Current, mA/um^2  ']);
model.result(pg).name('DriftCurrent');

model.result(pg).feature('surf1').set('rangecoloractive', 'on');
model.result(pg).feature('surf1').set('rangecolormin', '0');
model.result(pg).feature('surf1').set('rangecolormax', '550');

model.result(pg).feature.create('arws1', 'ArrowSurface');
model.result(pg).feature('arws1').setIndex('expr', '-ISpinDriftx', 0);
model.result(pg).feature('arws1').setIndex('expr', '-ISpinDrifty', 1);
model.result(pg).feature('arws1').set('arrowlength', 'logarithmic');
end



model.result(pg).run;
h(j)=figure(j); mphplot(model,pg,'rangenum',1);
axis([-2.5,2.5,-2.5,2.5])
xlabel('x-coordinate, \mum','FontSize', 18)
ylabel('y-coordinate, \mum','FontSize', 18)
title({['Drift Spin Current, , uA/um^2']; commonTitle}, 'FontSize',16)
F(j11,j)=getframe(h(j));
 MyTitle{j}='SpinDrift';

end

if SaveToFile

    for j2=1:j
        F1(:,1)=F(:,j2);
    save([myPath MyFile MyTitle{j2} 'F1' ],'F1');  %save images frame
  clear F1;
    end
end

close all;
h1=figure(1);movie(h1,F(:,1),100,3);

disp('All OK');model
return