% COMSOL Multiphysics Model M-file
% Generated by COMSOL 3.4 (COMSOL 3.4.0.250, $Date: 2008/01/24 15:34:30 $)

clear

%Loop Parameters
%tmin = 1;
%tmax = 2;
smin = 1;
smax = 15;

%Defining constants
T = 300;
kb = 1.38E-23;
e = 1.602E-19;
epsilon0 = 8.85E-12;
epsilon1 = 11.7*epsilon0;
epsilon2 = 4*epsilon0;
epsilon3 = 4*epsilon0;
epsilon4 = 78*epsilon0;

%Physical parameters
L = 10E-6;     %Length of Wire
r = 10E-9;      %Radius of wire
d = 1E-9;       %Thickness of Oxide
f = 1E-9;      %Thickness of functional layer
liq = 400E-9;   %Thickness of liquid

%Doping parameters
n0 = 1.2E21;    
ni = 1E16;
p0 = ni.^2./n0;

%mobility
un = 1300/1E4;
up = 500/1E4;

%Surface Charge
Sigma12 = 0;
Sigma23 = 0;

Q1 = 0.05E14;
pKa1 = 7.21;
pH = [7, linspace(4,8.5,smax)];
Sigma34 = [Q1 ./ (1 + 10.^(pH - pKa1))];

%Charge Density
rho3 = 0;

%Ion concentration
c0 = [10*6.022E23];
Z = 1;

for j = smin:smax+1
          
flclear fem

% COMSOL version
clear vrsn
vrsn.name = 'COMSOL 3.4';
vrsn.ext = '';
vrsn.major = 0;
vrsn.build = 250;
vrsn.rcs = '$Name:  $';
vrsn.date = '$Date: 2008/01/24 15:34:30 $';
fem.version = vrsn;

% Geometry
g1=circ2(r,'base','center','pos',{'0','0'},'rot','0');
g2=circ2(r+d,'base','center','pos',{'0','0'},'rot','0');
g3=circ2(r+d+f,'base','center','pos',{'0','0'},'rot','0');
g4=circ2(r+d+f+liq,'base','center','pos',{'0','0'},'rot','0');

% Constants
fem.const = {'epsilon0',epsilon0, ...
  'epsilon1',epsilon1, ...
  'epsilon2',epsilon2, ...
  'epsilon3',epsilon3, ...
  'epsilon4',epsilon4, ...
  'e',e, ...
  'n0',n0, ...
  'L', L, ...
  'un', un, ...
  'up', up, ...
  'ni',ni, ...
  'p0',p0, ...
  'kb',kb, ...
  'T',T, ...
  'Z',Z, ...
  'pi',pi, ...
  'c0',c0, ...
  'rho3',rho3, ...
  'sigma12',Sigma12, ...
  'sigma23',Sigma23, ...
  'sigma34',Sigma34(j).*e.*1E2.*1E2};

% Geometry

% Analyzed geometry
clear s
s.objs={g1,g2,g3,g4};
s.name={'Wire','Oxide','Func','Liquid'};
s.tags={'g1','g2','g3','g4'};

fem.draw=struct('s',s);
fem.geom=geomcsg(fem);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',5);

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'Poisson';
appl.border = 'on';
appl.assignsuffix = '_poeq';
clear bnd
bnd.g = {0,'-sigma34','-sigma23','-sigma12'};
bnd.type = {'dir','neu','neu','neu'};
bnd.ind = [1,1,2,2,3,3,4,4,1,2,3,4,4,3,2,1];
appl.bnd = bnd;
clear equ
equ.f = {'-2*Z*e*c0*sinh(Z*e*u/(kb*T))','rho3',0,'-e*(n0-p0-n0*exp(-e*u/(kb*T))+p0*exp(e*u/(kb*T)))'};
equ.c = {'epsilon4','epsilon3','epsilon2','epsilon1'};
equ.ind = [1,2,3,4];
appl.equ = equ;
fem.appl{1} = appl;
fem.frame = {'ref'};
fem.border = 1;
clear units;
units.basesystem = 'SI';
fem.units = units;

% Descriptions
clear descr
descr.const= {'n0','Electron density','e','Electron Charge','p0','Hole density','sigma23','oxide-func','epsilon4','Dielectric constant of liquid','T','Temperature','epsilon3','Dielectric constant of func','sigma12','Surface charge silicon-oxide','epsilon2','Dielectric constant of oxide','rho3','Charge distribution of func','sigma34','func-liquid','c0','Ion concentration','kb','Boltzmann constant','ni','Intrinsic charge','epsilon1','Dielectric constant of silicon','epsilon0','Dieelctric constant','Z','Electronic charge'};
fem.descr = descr;

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

% Solve problem
fem.sol=femstatic(fem, ...
                  'solcomp',{'u'}, ...
                  'outcomp',{'u'});

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

G(j)=postint(fem,'1/L*e*(un*n0*exp(-e*u/(kb*T))+up*p0*exp(e*u/(kb*T)))', ...
           'unit','', ...
           'dl',[4]);
       

       
       PercentageCompleted = ((j)/(smax))*100
end

for j=2:smax+1
    GG(j-1) = (G(j)-G(1))/G(1)*100;
    pH2(j-1) = pH(j)
end


plot(pH2,GG,'x-b')
xlabel('pH')
ylabel('\DeltaR / R_0 [%]')
title('Change in resistance as a function of pH')
legend('Theoretical Model')