neumannEquationBC.new.h

Go to the documentation of this file.

c File generated by ogmg/neumann.maple
c Fourth-order BC for mixed boundary condition using the normal derivative of the PDE 
c BC:    an1*ur+an2*us+an0*u = g 
c PDE:   c11*urr + c12*urs + c22*uss + c1*ur + c2*us + c0 *u = f 
c PDE.r: urrr + br2*urr = b0*u + b1*us + b3*usss + bf  

c real alpha1,alpha2,a1,a2,a0
c real rxi,ryi,sxi,syi,rxr,rxs,sxr,sxs,ryr,rys,syr,sys
c real rxx,ryy,sxx,syy
c real rxrr,rxrs,rxss,ryrr,ryrs,ryss
c real sxrr,sxrs,sxss,syrr,syrs,syss
c real rxx,ryy,sxx,syy
c real rxxr,ryyr,rxxs,ryys, sxxr,syyr,sxxs,syys
c real rxNormI,rxNormIs,rxNormIss,rxNormIr,rxNormIrr
c real n1,n1s,n1ss,n2,n2s,n2ss
c real an1,an1s,an1ss,an2,an2s,an2ss 
c real an1r,an1rr,an2r,an2rs 
c real ff,ffs,ffr,g,gs,gss 
c real c11,c11r,c11s,c12,c12r,c12s,c22,c22r,c22s,c1,c1r,c1s,c2,c2r,c2s,c0,c0r,c0s
c real b0,b1,b2,b3,bf,br2

#beginMacro fourthOrderNeumannEquationBCNew(DIR,DIM)
 #If #DIM == "3" 
   stop 45
 #End
 rxi=rx(i1,i2,i3)
 ryi=ry(i1,i2,i3)
 sxi=sx(i1,i2,i3)
 syi=sy(i1,i2,i3)

 rxr=rxr2(i1,i2,i3)
 rxs=rxs2(i1,i2,i3)
 ryr=ryr2(i1,i2,i3)
 rys=rys2(i1,i2,i3)

 sxr=sxr2(i1,i2,i3)
 sxs=sxs2(i1,i2,i3)
 syr=syr2(i1,i2,i3)
 sys=sys2(i1,i2,i3)

 rxx=rxx22(i1,i2,i3)
 ryy=ryy22(i1,i2,i3)
 rxrr=rxrr2(i1,i2,i3)
 rxrs=rxrs2(i1,i2,i3)
 rxss=rxss2(i1,i2,i3)
 ryrr=ryrr2(i1,i2,i3)
 ryrs=ryrs2(i1,i2,i3)
 ryss=ryss2(i1,i2,i3)

 sxx=sxx22(i1,i2,i3)
 syy=syy22(i1,i2,i3)
 sxrr=sxrr2(i1,i2,i3)
 sxrs=sxrs2(i1,i2,i3)
 sxss=sxss2(i1,i2,i3)
 syrr=syrr2(i1,i2,i3)
 syrs=syrs2(i1,i2,i3)
 syss=syss2(i1,i2,i3)

 rxxr=rxi*rxrr+rxr*rxr + sxi*rxrs + sxr*rxs
 ryyr=ryi*ryrr+ryr*ryr + syi*ryrs + syr*rys

 rxxs=rxi*rxrs+rxs*rxr + sxi*rxss + sxs*rxs
 ryys=ryi*ryrs+rys*ryr + syi*ryss + sys*rys

 sxxr=rxi*sxrr+rxr*sxr + sxi*sxrs + sxr*sxs
 syyr=ryi*syrr+ryr*syr + syi*syrs + syr*sys

 sxxs=rxi*sxrs+rxs*sxr + sxi*sxss + sxs*sxs
 syys=ryi*syrs+rys*syr + syi*syss + sys*sys

 alpha1=a1*nsign  ! nsign=2*side-1
 alpha2=a1*nsign

 #If #DIR == "R"
  rxNormI=1./sqrt(rxi**2+ryi**2)
  rxNormIs=-(rxi*rxs+ryi*rys)*rxNormI**3
  rxNormIss=-(rxi*rxss+ryi*ryss+rxs*rxs+rys*rys)*rxNormI**3 -3.*(rxi*rxs+ryi*rys)*rxNormI**2*rxNormIs

  n1=rxi*rxNormI
  n1s=rxs*rxNormI + rxi*rxNormIs
  n1ss=rxss*rxNormI + 2.*rxs*rxNormIs + rxi*rxNormIss
  n2=ryi*rxNormI
  n2s=rys*rxNormI + ryi*rxNormIs
  n2ss=ryss*rxNormI + 2.*rys*rxNormIs + ryi*rxNormIss
  an1=alpha1*(n1*rxi+n2*ryi)
  an2=alpha2*(n1*sxi+n2*syi)
  an1s=alpha1*(n1*rxs+n2*rys+n1s*rxi+n2s*ryi)
  an1ss=alpha1*(n1*rxss+n2*ryss+2.*(n1s*rxs+n2s*rys)+n1ss*rxi+n2ss*ryi)
  an2s=alpha2*(n1*sxs+n2*sys + n1s*sxi+n2s*syi)
  an2ss=alpha2*(n1*sxss+n2*syss + 2.*(n1s*sxs+n2s*sys) + n1ss*sxi+n2ss*syi)
 #Else
  sxNormI=1./sqrt(sxi**2+syi**2)
  sxNormIr=-(sxi*sxr+syi*syr)*sxNormI**3
  sxNormIrr=-(sxi*sxrr+syi*syrr+sxr*sxr+syr*syr)*sxNormI**3 -3.*(sxi*sxr+syi*syr)*sxNormI**2*sxNormIr

  n1=sxi*sxNormI
  n1r=sxr*sxNormI + sxi*sxNormIr
  n1rr=sxrr*sxNormI + 2.*sxr*sxNormIr + sxi*sxNormIrr
  n2=syi*sxNormI
  n2r=syr*sxNormI + syi*sxNormIr
  n2rr=syrr*sxNormI + 2.*syr*sxNormIr + syi*sxNormIrr
  an1=alpha1*(n1*rxi+n2*ryi)
  an2=alpha2*(n1*sxi+n2*syi)
  an1r=alpha1*(n1*rxr+n1r*rxi + n2*ryr+n2r*ryi)
  an1rr=alpha1*(n1*rxrr+n2*ryrr+ 2.*(n1r*rxr+n2r*ryr) + n1rr*rxi+n2rr*ryi)
  an2r=alpha2*(n1*sxr+n1r*sxi + n2*syr+n2r*syi)
  an2rr=alpha2*(n1*sxrr+n2*syrr+ 2.*(n1r*sxr+n2r*syr) + n1rr*sxi+n2rr*syi)
 #End

 c11=rxi**2+ryi**2
 c11r=2.*(rxi*rxr+ryi*ryr)
 c11s=2.*(rxi*rxs+ryi*rys)
 c12=2.*(rxi*sxi+ryi*syi)
 c12r=2.*(rxr*sxi+rxi*sxr + ryr*syi+ ryi*syr)
 c12s=2.*(rxs*sxi+rxi*sxs + rys*syi+ ryi*sys)
 c22=sxi**2+syi**2
 c22r=2.*(sxi*sxr+syi*syr)
 c22s=2.*(sxi*sxs+syi*sys)
 c1=rxx+ryy
 c1r=rxxr+ryyr
 c1s=rxxs+ryys
 c2=sxx+syy
 c2r=sxxr+syyr
 c2s=sxxs+syys
 c0=0.
 c0r=0.
 c0s=0.

c Coefficients of u(i1-2,i2,i3) u(i1-1,i2,i3) u(i1,i2,i3) u(i1+1,i2,i3) u(i1+2,i2,i3)
c um2=-1/2/dr**3
c um1=1/dr**3+1/dr**2*br2
c u0 =-b0-2/dr**2*br2
c up1=1/dr**2*br2-1/dr**3
c up2=1/2/dr**3
 #If #DIR == "R"
 b0=-(c12*c1*an1**2*an2*c0+c1s*an1*a0*c12**2*an2-c22s*an1*c12**2*an1s*a0-c0s*an1**2*c12**2*an2+2*an2s*c12**3*an1s*a0+c12**2*an1ss*a0*c22*an1-2*c12**2*an2s*c1*an1*a0-c0r*an1**3*c11*c22+c0r*an1**2*c11*c12*an2+c1r*an1**2*c11*a0*c22-c22r*an1**2*c11*c1*a0-c1r*an1*c11*a0*c12*an2-c12r*an1**2*c11*an2*c0+c12r*an1*c11*an2*c1*a0-c2*an1**2*c11*an2*c0+c2*an1*c11*an2*c1*a0+c22r*an1**3*c11*c0-c22**2*c11*an1ss*a0*an1+c22*c11*c1*an1*a0**2-c22*c11*c12*an1s*a0**2-2*c22*c11*an2s*c0*an1**2+2*c22**2*c11*an1s**2*a0-c12*c1s*an1**2*a0*c22-c12*c22s*an1**3*c0-c12*c2*an1**3*c0-c12*c1**2*an1*an2*a0+c12*c2*an1**2*c1*a0+2*c12**2*an2s*c0*an1**2-an1ss*a0*c12**3*an2-2*c12**2*an1s*an2*c0*an1+c12**3*an1s*a0**2+2*c22*c11*an2s*c1*an1*a0-2*c22*c11*an2s*c12*an1s*a0+c12*c12s*an1*an1s*a0*c22-c12*c12s*an1*an2*c1*a0+c12*c0s*an1**3*c22+c12*c1*an1*an1s*a0*c22+c12**2*a0*c0*an1**2-2*c12**2*an1s**2*a0*c22-c12**2*c1*an1*a0**2+2*c12**2*an1s*an2*c1*a0+c22r*an1*c11*c12*an1s*a0-c12r*an1*c11*an1s*a0*c22-c0*an1*c11*a0*c12*an2-c2*an1*c11*an1s*a0*c22+c22*c11*an1ss*a0*c12*an2-2*c22*c11*an1s*an2*c1*a0+2*c22*c11*an1s*an2*c0*an1+c12*c12s*an1**2*an2*c0+c12*c22s*an1**2*c1*a0-c2*an1*c12**2*an1s*a0)/an1**2/(c12*an2-c22*an1)/c11**2

 b1=1/an1**2*(2*c12**2*an2s*c1*an1*an2+c12*c12s*an1**2*an2s*c22+c12*c12s*an1**2*a0*c22+c22r*an1**2*c11*c12*a0-c22r*an1*c11*c12*an1s*an2+c22r*an1**2*c11*c1*an2+c1r*an1*c11*an2**2*c12-c1r*an1**2*c11*an2*c22-c12*c12s*an1**2*an2*c2-c12*c12s*an1*an1s*an2*c22-c2r*an1**2*c11*c12*an2+c12**2*a0*c1*an1*an2+3*c12**2*an1s*an2*c2*an1-2*c12**2*an1s*an2s*an1*c22+2*c12**3*an2s**2*an1-c12r*an1*c11*c1*an2**2-c2*an1*c11*c1*an2**2+c0*an1*c11*an2**2*c12-c0*an1**2*c11*an2*c22+c12r*an1**2*c11*an2*c2-3*c22*c11*an2s*c12*a0*an1+2*c22*c11*an2s*c12*an1s*an2+c2**2*an1**2*c11*an2-c22*c11*c12*a0**2*an1+2*c22*c11*an1s*c1*an2**2-c22*c11*an1ss*an2**2*c12-2*c22*c11*an2s*c1*an1*an2+c22*c11*an2ss*an1*c12*an2+c22*c11*a0*c12*an1s*an2-c22*c11*an1s*an2*c2*an1+2*c22**2*c11*an1s*an2s*an1+2*c22**2*c11*an1s*a0*an1+c22**2*c11*an1ss*an2*an1+3*an2s*c12**3*a0*an1-2*an2s*c12**3*an1s*an2-an2ss*an1*c12**3*an2+2*c12**2*an1s**2*an2*c22-2*c12**2*an1s*c1*an2**2+c12*c1**2*an1*an2**2-c12*c0*an1**3*c22+c12*c22s*an1**3*c2-2*c12**2*an1s*a0*an1*c22-c12**2*an1ss*an2*c22*an1-a0*c12**3*an1s*an2+c12**2*an2ss*an1**2*c22-2*c12**2*a0*c2*an1**2-3*c12**2*an2s*c2*an1**2-2*c12*c1*an1**2*an2*c2+c12*c1*an1**2*an2s*c22+c12*c1*an1**2*a0*c22-2*c22*c11*c12*an2s**2*an1+c22*c11*an2s*c2*an1**2-c12*c22s*an1**2*c1*an2+c12*c1s*an1**2*an2*c22-c22*c11*a0*c1*an1*an2-c12r*an1**2*c11*an2s*c22-c12r*an1**2*c11*a0*c22+c12r*an1*c11*an1s*an2*c22+c22r*an1**2*c11*c12*an2s+c12**3*a0**2*an1-c12*c2s*an1**3*c22-c22**2*c11*an2ss*an1**2-2*c22**2*c11*an1s**2*an2-c22r*an1**3*c11*c2+c2r*an1**3*c11*c22-c12*c1*an1*an1s*an2*c22+c12*c12s*an1*c1*an2**2+c22s*an1*c12**2*an1s*an2-c1s*an1*an2**2*c12**2+c2s*an1**2*c12**2*an2-c22s*an1**2*c12**2*a0+c0*an1**2*c12**2*an2-c22s*an1**2*c12**2*an2s+c12*c2**2*an1**3+an1ss*an2**2*c12**3)/(c12*an2-c22*an1)/c11**2

 br2=1/an1*(-c11s*an1*c12**2*an2+2*c12**2*an2s*c11*an1+c12**2*a0*c11*an1+2*c12*c1*an1*an2*c11+c12*c12s*an1*an2*c11-c12*c2*an1**2*c11-c12*c22s*an1**2*c11+c12*c11s*an1**2*c22-2*c22*c11**2*an2s*an1-c22*c11**2*a0*an1+2*c22*c11**2*an1s*an2-c1*an1**2*c11*c22-c2*an1*c11**2*an2+c22r*an1**2*c11**2-c12r*an1*c11**2*an2+c11r*an1*c11*c12*an2-c11r*an1**2*c11*c22-2*c12**2*an1s*an2*c11)/(c12*an2-c22*an1)/c11**2

 b3=-(an2*c12**2-an2*c22*c11-c22*an1*c12)/an1/c11**2

 bf=1/an1**2*(2*c12**2*an1s*an1*c22-c12*c1*an1**2*c22-c12*c12s*an1**2*c22+2*c22*c11*an2s*c12*an1+c22*c11*a0*c12*an1-2*c22**2*c11*an1s*an1+c2*an1**2*c11*c22+c12r*an1**2*c11*c22-c22r*an1**2*c11*c12-2*an2s*c12**3*an1-a0*c12**3*an1+c2*an1**2*c12**2+c22s*an1**2*c12**2)/(c12*an2-c22*an1)/c11**2*gs+1/an1**2*(-c22*c11*an1*c12*an2+an1*c12**3*an2-c12**2*an1**2*c22+c22**2*c11*an1**2)/(c12*an2-c22*an1)/c11**2*gss+1/an1**2*(-2*c12**2*an1s*an2*an1+c12*c1*an1**2*an2+c12*c12s*an1**2*an2+2*c22*c11*an1s*an2*an1-c22*c11*a0*an1**2-2*c22*c11*an2s*an1**2-c2*an1**2*c11*an2+c22r*an1**3*c11-c12r*an1**2*c11*an2+c12**2*a0*an1**2+2*c12**2*an2s*an1**2-c12*c2*an1**3-c12*c22s*an1**3)/(c12*an2-c22*an1)/c11**2*ff+1/an1**2*(-an1**2*c12**2*an2+c12*an1**3*c22)/(c12*an2-c22*an1)/c11**2*ffs+1/an1**2*(-an1**3*c11*c22+an1**2*c11*c12*an2)/(c12*an2-c22*an1)/c11**2*ffr+1/an1**2*(-2*c12**2*an2s*c1*an1-c22s*an1*c12**2*an1s-c2*an1*c12**2*an1s+c1s*an1*c12**2*an2-an1ss*c12**3*an2-2*c12**2*an1s**2*c22+c22r*an1*c11*c12*an1s-c12r*an1*c11*an1s*c22+c12r*an1*c11*an2*c1-c22r*an1**2*c11*c1-c0*an1*c11*c12*an2+2*an2s*c12**3*an1s-c2*an1*c11*an1s*c22-c12*c1s*an1**2*c22+c12*c22s*an1**2*c1+c2*an1*c11*an2*c1+c12**2*an1ss*c22*an1-c12**2*a0*c1*an1+c22*c11*an1ss*c12*an2-2*c22*c11*an1s*an2*c1+c22*c11*a0*c1*an1-c22**2*c11*an1ss*an1+c0*an1**2*c11*c22+c12*c12s*an1*an1s*c22-c12*c1**2*an1*an2-c12*c12s*an1*an2*c1+c12*c2*an1**2*c1-2*c22*c11*an2s*c12*an1s+a0*c12**3*an1s+c12*c1*an1*an1s*c22+2*c12**2*an1s*an2*c1+2*c22**2*c11*an1s**2+2*c22*c11*an2s*c1*an1-c22*c11*a0*c12*an1s-c1r*an1*c11*c12*an2+c1r*an1**2*c11*c22)/(c12*an2-c22*an1)/c11**2*g

 #Else
 b0=-(c12*c2*an2**2*an1*c0+c2r*an2*a0*c12**2*an1-c11r*an2*c12**2*an2r*a0-c0r*an2**2*c12**2*an1+2*an1r*c12**3*an2r*a0+c12**2*an2rr*a0*c11*an2-2*c12**2*an1r*c2*an2*a0-c0s*an2**3*c22*c11+c0s*an2**2*c22*c12*an1+c2s*an2**2*c22*a0*c11-c11s*an2**2*c22*c2*a0-c2s*an2*c22*a0*c12*an1-c12s*an2**2*c22*an1*c0+c12s*an2*c22*an1*c2*a0-c1*an2**2*c22*an1*c0+c1*an2*c22*an1*c2*a0+c11s*an2**3*c22*c0-c11**2*c22*an2rr*a0*an2+c11*c22*c2*an2*a0**2-c11*c22*c12*an2r*a0**2-2*c11*c22*an1r*c0*an2**2+2*c11**2*c22*an2r**2*a0-c12*c2r*an2**2*a0*c11-c12*c11r*an2**3*c0-c12*c1*an2**3*c0-c12*c2**2*an2*an1*a0+c12*c1*an2**2*c2*a0+2*c12**2*an1r*c0*an2**2-an2rr*a0*c12**3*an1-2*c12**2*an2r*an1*c0*an2+c12**3*an2r*a0**2+2*c11*c22*an1r*c2*an2*a0-2*c11*c22*an1r*c12*an2r*a0+c12*c12r*an2*an2r*a0*c11-c12*c12r*an2*an1*c2*a0+c12*c0r*an2**3*c11+c12*c2*an2*an2r*a0*c11+c12**2*a0*c0*an2**2-2*c12**2*an2r**2*a0*c11-c12**2*c2*an2*a0**2+2*c12**2*an2r*an1*c2*a0+c11s*an2*c22*c12*an2r*a0-c12s*an2*c22*an2r*a0*c11-c0*an2*c22*a0*c12*an1-c1*an2*c22*an2r*a0*c11+c11*c22*an2rr*a0*c12*an1-2*c11*c22*an2r*an1*c2*a0+2*c11*c22*an2r*an1*c0*an2+c12*c12r*an2**2*an1*c0+c12*c11r*an2**2*c2*a0-c1*an2*c12**2*an2r*a0)/an2**2/(c12*an1-c11*an2)/c22**2

 b1=1/an2**2*(2*c12**2*an1r*c2*an2*an1+c12*c12r*an2**2*an1r*c11+c12*c12r*an2**2*a0*c11+c11s*an2**2*c22*c12*a0-c11s*an2*c22*c12*an2r*an1+c11s*an2**2*c22*c2*an1+c2s*an2*c22*an1**2*c12-c2s*an2**2*c22*an1*c11-c12*c12r*an2**2*an1*c1-c12*c12r*an2*an2r*an1*c11-c1s*an2**2*c22*c12*an1+c12**2*a0*c2*an2*an1+3*c12**2*an2r*an1*c1*an2-2*c12**2*an2r*an1r*an2*c11+2*c12**3*an1r**2*an2-c12s*an2*c22*c2*an1**2-c1*an2*c22*c2*an1**2+c0*an2*c22*an1**2*c12-c0*an2**2*c22*an1*c11+c12s*an2**2*c22*an1*c1-3*c11*c22*an1r*c12*a0*an2+2*c11*c22*an1r*c12*an2r*an1+c1**2*an2**2*c22*an1-c11*c22*c12*a0**2*an2+2*c11*c22*an2r*c2*an1**2-c11*c22*an2rr*an1**2*c12-2*c11*c22*an1r*c2*an2*an1+c11*c22*an1rr*an2*c12*an1+c11*c22*a0*c12*an2r*an1-c11*c22*an2r*an1*c1*an2+2*c11**2*c22*an2r*an1r*an2+2*c11**2*c22*an2r*a0*an2+c11**2*c22*an2rr*an1*an2+3*an1r*c12**3*a0*an2-2*an1r*c12**3*an2r*an1-an1rr*an2*c12**3*an1+2*c12**2*an2r**2*an1*c11-2*c12**2*an2r*c2*an1**2+c12*c2**2*an2*an1**2-c12*c0*an2**3*c11+c12*c11r*an2**3*c1-2*c12**2*an2r*a0*an2*c11-c12**2*an2rr*an1*c11*an2-a0*c12**3*an2r*an1+c12**2*an1rr*an2**2*c11-2*c12**2*a0*c1*an2**2-3*c12**2*an1r*c1*an2**2-2*c12*c2*an2**2*an1*c1+c12*c2*an2**2*an1r*c11+c12*c2*an2**2*a0*c11-2*c11*c22*c12*an1r**2*an2+c11*c22*an1r*c1*an2**2-c12*c11r*an2**2*c2*an1+c12*c2r*an2**2*an1*c11-c11*c22*a0*c2*an2*an1-c12s*an2**2*c22*an1r*c11-c12s*an2**2*c22*a0*c11+c12s*an2*c22*an2r*an1*c11+c11s*an2**2*c22*c12*an1r+c12**3*a0**2*an2-c12*c1r*an2**3*c11-c11**2*c22*an1rr*an2**2-2*c11**2*c22*an2r**2*an1-c11s*an2**3*c22*c1+c1s*an2**3*c22*c11-c12*c2*an2*an2r*an1*c11+c12*c12r*an2*c2*an1**2+c11r*an2*c12**2*an2r*an1-c2r*an2*an1**2*c12**2+c1r*an2**2*c12**2*an1-c11r*an2**2*c12**2*a0+c0*an2**2*c12**2*an1-c11r*an2**2*c12**2*an1r+c12*c1**2*an2**3+an2rr*an1**2*c12**3)/(c12*an1-c11*an2)/c22**2

 br2=1/an2*(-c22r*an2*c12**2*an1+2*c12**2*an1r*c22*an2+c12**2*a0*c22*an2+2*c12*c2*an2*an1*c22+c12*c12r*an2*an1*c22-c12*c1*an2**2*c22-c12*c11r*an2**2*c22+c12*c22r*an2**2*c11-2*c11*c22**2*an1r*an2-c11*c22**2*a0*an2+2*c11*c22**2*an2r*an1-c2*an2**2*c22*c11-c1*an2*c22**2*an1+c11s*an2**2*c22**2-c12s*an2*c22**2*an1+c22s*an2*c22*c12*an1-c22s*an2**2*c22*c11-2*c12**2*an2r*an1*c22)/(c12*an1-c11*an2)/c22**2

 b3=-(an1*c12**2-an1*c11*c22-c11*an2*c12)/an2/c22**2

 bf=1/an2**2*(2*c12**2*an2r*an2*c11-c12*c2*an2**2*c11-c12*c12r*an2**2*c11+2*c11*c22*an1r*c12*an2+c11*c22*a0*c12*an2-2*c11**2*c22*an2r*an2+c1*an2**2*c22*c11+c12s*an2**2*c22*c11-c11s*an2**2*c22*c12-2*an1r*c12**3*an2-a0*c12**3*an2+c1*an2**2*c12**2+c11r*an2**2*c12**2)/(c12*an1-c11*an2)/c22**2*gr+1/an2**2*(-c11*c22*an2*c12*an1+an2*c12**3*an1-c12**2*an2**2*c11+c11**2*c22*an2**2)/(c12*an1-c11*an2)/c22**2*grr+1/an2**2*(-2*c12**2*an2r*an1*an2+c12*c2*an2**2*an1+c12*c12r*an2**2*an1+2*c11*c22*an2r*an1*an2-c11*c22*a0*an2**2-2*c11*c22*an1r*an2**2-c1*an2**2*c22*an1+c11s*an2**3*c22-c12s*an2**2*c22*an1+c12**2*a0*an2**2+2*c12**2*an1r*an2**2-c12*c1*an2**3-c12*c11r*an2**3)/(c12*an1-c11*an2)/c22**2*ff+1/an2**2*(-an2**2*c12**2*an1+c12*an2**3*c11)/(c12*an1-c11*an2)/c22**2*ffr+1/an2**2*(-an2**3*c22*c11+an2**2*c22*c12*an1)/(c12*an1-c11*an2)/c22**2*ffs+1/an2**2*(-2*c12**2*an1r*c2*an2-c11r*an2*c12**2*an2r-c1*an2*c12**2*an2r+c2r*an2*c12**2*an1-an2rr*c12**3*an1-2*c12**2*an2r**2*c11+c11s*an2*c22*c12*an2r-c12s*an2*c22*an2r*c11+c12s*an2*c22*an1*c2-c11s*an2**2*c22*c2-c0*an2*c22*c12*an1+2*an1r*c12**3*an2r-c1*an2*c22*an2r*c11-c12*c2r*an2**2*c11+c12*c11r*an2**2*c2+c1*an2*c22*an1*c2+c12**2*an2rr*c11*an2-c12**2*a0*c2*an2+c11*c22*an2rr*c12*an1-2*c11*c22*an2r*an1*c2+c11*c22*a0*c2*an2-c11**2*c22*an2rr*an2+c0*an2**2*c22*c11+c12*c12r*an2*an2r*c11-c12*c2**2*an2*an1-c12*c12r*an2*an1*c2+c12*c1*an2**2*c2-2*c11*c22*an1r*c12*an2r+a0*c12**3*an2r+c12*c2*an2*an2r*c11+2*c12**2*an2r*an1*c2+2*c11**2*c22*an2r**2+2*c11*c22*an1r*c2*an2-c11*c22*a0*c12*an2r-c2s*an2*c22*c12*an1+c2s*an2**2*c22*c11)/(c12*an1-c11*an2)/c22**2*g

 #End
#endMacro
c Solution for ur4=ga, urrr2+br2*urr2 = gb is 
c u(i1-1,i2,i3) = 2*br2*drn/(-3+br2*drn)*u(i1,i2,i3)+(-3-br2*drn)/(-3+br2*drn)*u(i1+1,i2,i3)+(6*ga*drn+gb*drn**3)/(-3+br2*drn)
c u(i1-2,i2,i3) = u(i1+2,i2,i3)+16*br2*drn/(-3+br2*drn)*u(i1,i2,i3)-16*br2*drn/(-3+br2*drn)*u(i1+1,i2,i3)+(12*ga*drn+8*gb*drn**3+12*ga*drn**2*br2)/(-3+br2*drn)

c Solution for ur4=ga, urrr2+br2*urr2 = gb is 
c u(i1+1,i2,i3) = 2*br2*drn/(3+br2*drn)*u(i1,i2,i3)+(3-br2*drn)/(3+br2*drn)*u(i1-1,i2,i3)+(6*ga*drn+gb*drn**3)/(3+br2*drn)
c u(i1+2,i2,i3) = u(i1-2,i2,i3)+16*br2*drn/(3+br2*drn)*u(i1,i2,i3)-16*br2*drn/(3+br2*drn)*u(i1-1,i2,i3)+(12*ga*drn+8*gb*drn**3-12*ga*drn**2*br2)/(3+br2*drn)

c Solution for ur4=ga, urrr2+br2*urr2 = gb, ca*u(i1-2)+cb*u(i1-1)+cc*u = gc is 
c u(i1-2,i2,i3) = (-3*cc+cc*br2*drn+2*br2*drn*cb)/(-3*cc+cc*br2*drn+16*br2*drn*ca+2*br2*drn*cb)*u(i1+2,i2,i3)+(-16*cc*br2*drn-16*br2*drn*cb)/(-3*cc+cc*br2*drn+16*br2*drn*ca+2*br2*drn*cb)*u(i1+1,i2,i3)+(16*br2*drn*gc+12*ga*drn**2*cc*br2+8*cc*gb*drn**3+12*cc*ga*drn+24*ga*drn**2*br2*cb)/(-3*cc+cc*br2*drn+16*br2*drn*ca+2*br2*drn*cb)
c u(i1,i2,i3) = -(br2*drn*ca-3*ca)/(-3*cc+cc*br2*drn+16*br2*drn*ca+2*br2*drn*cb)*u(i1+2,i2,i3)-(-16*br2*drn*ca-br2*drn*cb-3*cb)/(-3*cc+cc*br2*drn+16*br2*drn*ca+2*br2*drn*cb)*u(i1+1,i2,i3)-(3*gc+6*cb*ga*drn+12*br2*drn**2*ca*ga-br2*drn*gc+cb*gb*drn**3+8*ca*gb*drn**3+12*ca*ga*drn)/(-3*cc+cc*br2*drn+16*br2*drn*ca+2*br2*drn*cb)
c u(i1-1,i2,i3) = -2/(-3*cc+cc*br2*drn+16*br2*drn*ca+2*br2*drn*cb)*br2*drn*ca*u(i1+2,i2,i3)-1/(-3*cc+cc*br2*drn+16*br2*drn*ca+2*br2*drn*cb)*(3*cc+cc*br2*drn-16*br2*drn*ca)*u(i1+1,i2,i3)-1/(-3*cc+cc*br2*drn+16*br2*drn*ca+2*br2*drn*cb)*(-cc*gb*drn**3-6*cc*ga*drn+24*br2*drn**2*ca*ga-2*br2*drn*gc)

c Solution for ur4=ga, urrr2+br2*urr2 = gb, ca*u(i1-2)+cb*u(i1-1)+cc*u = gc is 
c u(i1+2,i2,i3) = -(16*br2*drn*ca-3*cc-cc*br2*drn)/cc/(3+br2*drn)*u(i1-2,i2,i3)-(16*cc*br2*drn+16*br2*drn*cb)/cc/(3+br2*drn)*u(i1-1,i2,i3)-(-12*cc*ga*drn+12*ga*drn**2*cc*br2-16*br2*drn*gc-8*cc*gb*drn**3)/cc/(3+br2*drn)
c u(i1,i2,i3) = -ca/cc*u(i1-2,i2,i3)-cb/cc*u(i1-1,i2,i3)+gc/cc
c u(i1+1,i2,i3) = -2*br2*drn*ca/cc/(3+br2*drn)*u(i1-2,i2,i3)-(-3*cc+cc*br2*drn+2*br2*drn*cb)/cc/(3+br2*drn)*u(i1-1,i2,i3)-(-2*br2*drn*gc-6*cc*ga*drn-cc*gb*drn**3)/cc/(3+br2*drn)