{Fpe2dThDNA.PDE  - FPE PDESolns - RGN's advice f: FIRSTPARTS }
title 'Fpe2DThDNA: 2D FPE Thesis - DNA' {D.Bakewell 23:30 We/21/11/2001}
{Thesis ch 7 alterations for DNA polarisability etc}
SELECT
errlim=8e-5  regrid=on    firstparts
COORDINATES
cartesian (x, y)
VARIABLES 
p     {solute probability or concentration}
DEFINITIONS
a=0.1 b={200.1} 60.1   {rectangle horiz and vert dimnsions, respctivly}
d=10  d12=d/2 d15=d/5 d10=d/10 dA=5.05/2 d2=2*d  Vo=0.45
xc=d/2 xcg=d-dA  xl=xcg-dA xu=xcg+dA { xcg is ca. 1/2 gap frm + el rght eg}
 alm= {0.2} {1.9}  7.9 {8.723} {DNA polarisability - see Book III, p.90-1}
 r = 0.200 {0.020} {0.002}
Kdiff=0.2424/r  { = 2.245 = KbT*(10^12)/zeta }
Kfpe=18.76*alm*(1/r)*(Vo^2)/(d^3) {nrmlizd valu - see Book III, p.91}
Kss=6.159*alm*(Vo^2)/(d^2) {nrmlizd valu - see Book III, p.91}
 xh=(pi*x/d2)+(pi/4)  yh=pi*y/d2 { h = 'hat' } 
ThT=(ARCTAN(sin(xh)/sinh(yh))-ARCTAN(cos(xh)/sinh(yh)))
ThL=(LN(cosh(yh)+cos(xh))-LN(cosh(yh)-cos(xh))+
  LN(cosh(yh)+sin(xh))-LN(cosh(yh)-sin(xh)))
Dpos=cosh(2*yh)+cos(2*xh)  Dneg=cosh(2*yh)-cos(2*xh)
XiS=sinh(yh)*(cos(xh)/Dneg+sin(xh)/Dpos) {Xi sum}
XiD=cosh(yh)*(cos(xh)/Dpos-sin(xh)/Dneg)  {Xi diff}
GamX=2.*ThT*XiS+ThL*XiD GamY=2.*ThT*XiD-ThL*XiS
Vx=Kfpe*GamX  Vy=Kfpe*GamY
Areap=(2*dA)*1  {Areap=pi*(dA^2)/2}  pA=integral(p,2)/Areap
Gamdiff=dx(GamX) + dy(GamY)
GamdiffAv=integral(Gamdiff)/(d*(b-a))
GamdiffAv2=integral(Gamdiff,2)/Areap
dpAdt_p0=-Kfpe*GamdiffAv2
diffusion=Kdiff*div[grad(p)]
difAv=integral(diffusion,2)/Areap
{- steady state p(x,y) from analytical solution and flux monitoring -}
ThSq=4.*ThT^2+ThL^2        pss=EXP(Kss*ThSq) 
 Cint=integral(pss) 	{dble intgrl of p}
 IQx=integral(Vx) IQy=integral(Vy) IThSq=integral(ThSq)
ps=pss/Cint  psA=integral(ps,2)/Areap  {pInit=initial p} 
pInit=1/(d)*1/(b-a)  pIA=integral(pInit,2)/Areap 
pAtau=(psA-pIA)*(1.-exp(-1.))+pIA 
DpARel=(psA-pIA)/pIA
Jx= -(Kdiff*dx(p)-p*Vx)   {for monitoring Flux in X direction} 
Jy= -(Kdiff*dy(p)-p*Vy)    {for monitoring Flux in Y direction}
J= vector(p*Vx-Kdiff*dx(p),p*Vy-Kdiff*dy(p)) 
INITIAL VALUES
p=1/(d)*1/(b-a) 
EQUATIONS
Kdiff*div[grad(p)]-dx[p*Vx]-dy[Vy*p]=dt(p)      { 2D FPE equation }
BOUNDARIES
Region 1   natural(p) = 0  {walking outer boundary}
 start (xc-d12,a) line to (xc+d12,a) to (xc+d12,b) to (xc-d12,b) to finish 
 {start (xc-d,a) line to (xc+d,a) to (xc+d,b) to (xc-d,b) to finish}
Region 2 start  'Aint1' (xcg,a) line to (xcg-dA,a) to (xcg-dA,a+1)
 to (xcg+dA,a+1) to (xcg+dA,a) to finish 
{arc (center=xc,a)  to (xc+dA,a) line to finish}
FEATURE
start (xc-.1,a) line to (xc-.09,a) to (xc-.08,a) to (xc-.07,a) to (xc-.06,a)
 to (xc-.05,a) line to (xc-.04,a) to (xc-.03,a) to (xc-.02,a) to (xc-.01,a)
 to (xc,a) to (xc+.01,a) to (xc+.02,a) to (xc+.03,a) to (xc+.04,a) to
 (xc+.05,a) to (xc+.06,a) to (xc+.07,a) to (xc+.08,a) to (xc+.09,a) to (xc+.1,a)
TIME 0 to 2000 by 0.00000001     { time range and initial timestep (secs)}
MONITORS for cycle=1  contour (p) contour (LOG10(p)) 
{elevation (p) from (xc,a) to (xc,b)  as "X-section"}
elevation (LOG10(p)) from (xc,a) to (xc,b)  as "X-section"
elevation (LOG10(ps)) from (xc,a) to (xc,b)  as "X-section" 
contour(GamX) contour(GamY)
! contour(Vx) contour(Vy) {contour(LOG10(k*E2/Cdiff)) }
 contour (ps )  contour ((p-ps)/ps) 
 contour(LOG10(ABS(Gamdiff))) report(GamdiffAv) report(GamdiffAv2) report(dpAdt_p0) report(difAv)
 elevation(Gamdiff) from (xc,a) to (xc,b)  as "X-section"
 elevation(LOG10(Gamdiff)) from (xc-d12,a) to (xc+d12,a)  as "Y-section"
contour(LOG10(ps))report(r) report(Vo) report(alm) report(pIA) report(psA) report(pAtau)
HISTORIES history (p) at (xc,a) report(r) report(Vo) report (alm) report(xl) report(xu)
 history (pA) report(r) report(Vo) report(alm) report(pIA) report(psA) report(pAtau) 
 report(DpARel)
PLOTS     { plot prints a tag-delimited data list to the file "filename.TXT":}
 grid(x,y)
 history(pA) export format "#t#b#" file="c:\windows\desktop\0.45V60h0.200r7.9A2000s.txt"
end