diff options
| author | ryoji <ryoji@535fb057-194f-0410-b5a5-c63992f15602> | 1999-12-01 13:34:57 +0000 |
|---|---|---|
| committer | ryoji <ryoji@535fb057-194f-0410-b5a5-c63992f15602> | 1999-12-01 13:34:57 +0000 |
| commit | 5d27c1507b8402c8be8ec6433bbc5483db0203ec (patch) | |
| tree | ccef534f16a8818747f1faafb1d946e9fad79d45 /src | |
| parent | 4f2663df48d28c15194cde1de8e03a279534fb0b (diff) | |
This commit was generated by cvs2svn to compensate for changes in r2,
which included commits to RCS files with non-trunk default branches.
git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/RotatingDBHIVP/trunk@3 535fb057-194f-0410-b5a5-c63992f15602
Diffstat (limited to 'src')
| -rw-r--r-- | src/RotatingDBHIVP.F | 798 | ||||
| -rw-r--r-- | src/Stab3d.F | 277 | ||||
| -rw-r--r-- | src/cgauss.x | 57 | ||||
| -rw-r--r-- | src/gauss.x | 57 | ||||
| -rw-r--r-- | src/gij.x | 24 | ||||
| -rw-r--r-- | src/kij_axi.x | 62 | ||||
| -rw-r--r-- | src/kij_odd.x | 153 | ||||
| -rw-r--r-- | src/ksq2d_odd.x | 23 | ||||
| -rw-r--r-- | src/ksq_axi.x | 111 | ||||
| -rw-r--r-- | src/ksq_odd.x | 403 | ||||
| -rw-r--r-- | src/make.code.defn | 9 | ||||
| -rw-r--r-- | src/psi_1st_deriv.x | 20 | ||||
| -rw-r--r-- | src/psi_2nd_deriv.x | 71 |
13 files changed, 2065 insertions, 0 deletions
diff --git a/src/RotatingDBHIVP.F b/src/RotatingDBHIVP.F new file mode 100644 index 0000000..4f2e6b1 --- /dev/null +++ b/src/RotatingDBHIVP.F @@ -0,0 +1,798 @@ +/*@@ +c @file RotatingDBH.F +c @date 8 June 1999 +c @author Ryoji Takahashi +c @desc +c +c @enddesc +c@@*/ + +#include "cctk.h" +#include "cctk_parameters.h" +#include "cctk_arguments.h" + +c/*@@ +c @routine RotatingDBHIVP +c @date +c @author +c @desc +c +c @enddesc +c @calls +c @calledby +c @history +c +c @endhistory +c@@*/ + +c Need include file from Einstein +#include "CactusEinstein/Einstein/src/Einstein.h" + + subroutine RotatingDBHIVP(CCTK_FARGUMENTS) + + implicit none + + DECLARE_CCTK_FARGUMENTS + DECLARE_CCTK_PARAMETERS + +c Perhaps this and others should go into cctk.h + integer CCTK_Equals + + real*8 deta,dq,dphi + real*8 psi3d,detapsi3d,dqpsi3d,dphipsi3d,detaetapsi3d,detaqpsi3d, + $ detaphipsi3d,dqqpsi3d,dqpsi3d,dqphipsi3d,dphiphipsi3d + real*8, allocatable :: ac(:,:,:),ae(:,:,:),aw(:,:,:),an(:,:,:), + $ as(:,:,:),aq(:,:,:),ab(:,:,:),rhs(:,:,:),Ksq(:,:,:), + $ psisph(:,:,:),psip(:,:,:),psipp(:,:,:),detapsisph(:,:,:), + $ dqpsisph(:,:,:),dphipsisph(:,:,:),detaetapsisph(:,:,:), + $ detaqpsisph(:,:,:),detaphipsisph(:,:,:),dqqpsisph(:,:,:), + $ dqphipsisph(:,:,:),dphiphipsisph(:,:,:), + $ psi3d(:,:,:),detapsi3d(:,:,:),dqpsi3d(:,:,:), + $ dphipsi3d(:,:,:),detaetapsi3d(:,:,:),detaqpsi3d(:,:,:), + $ detaphipsi3d(:,:,:),dqqpsi3d(:,:,:),dqphipsi3d(:,:,:), + $ dphiphipsi3d(:,:,:) + real*8, allocatable :: etagrd(:),qgrd(:),phigrd(:) + real*8, allocatable :: eta(:,:,:),abseta(:,:,:),sign_eta(:,:,:), + $ q(:,:,:),phi(:,:,:) + real*8, allocatable ::r_BL(:,:,:),Delta(:,:,:),Rho(:,:,:), + $ beta_phi(:,:,:) + real*8 o1,o2,o3,o4,o5,o6,o7,o8,o9,o10,o11,o12,o13,o14,o15,o16,o17, + $ o18,o19,o20,o21,o22,o23,o24,o25,o26,o27,o28,o29,o30,o31,o32, + $ o33,o34,o35,o36,o37,o38,o39,o40,o41,o42,o43,o44,o45,o46,o47, + $ o48,o49,o50,o51,o52,o53,o54,o55,o56,o57,o58,o59,o60,o61,o62, + $ o63,o64,o65,o66,o67,o68,o69,o70,o71,o72,o73,o74,o75,o76,o77, + $ o78,o79,o80,o81,o82,o83,o84,o85,o86,o87,o88,o89,o90,o91,o92, + $ o93,o94,o95,o96,o97,o98,o99,o100,o101,o102,o103,o104,o105, + $ o106,o107,o108,o109,o110,o111,o112,o113,o114,o115,o116,o117, + $ o118,o119,o120,o121,o122,o123,o124 + real*8 t1,t2,t3,t4,t5,t6,t7,t8,t9,t10, + $ t11,t12,t13,t14,t15,t16 + real*8 gtil,dngtil,dnngtil,dnnngtil,dnnnngtil,dnnnnngtil, + $ dnnnnnngtil,dnnnnnnngtil + real*8 eta,abseta,q,phi,rhsmax,rmax,get3d,adm + real*8,parameter :: rbh_tol = 1.0d-7,rbh_eps = 1.0d-10 + integer,parameter :: itmax = 30 + real*8 r,pi + integer :: nx,ny,nz + integer i,j,k,ier,nquads,nocts,order,ntries + integer npoints,handle,ierror + integer i,j,k,it,ier,nquads,nocts,order + + conformal_state = CONFORMAL_METRIC + + pi = 4.0d0*atan(1.0d0) + +c Set up the grid spacings + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + +c Distorted Rotating BH parameters (for Extrinsic curveture!) +c + print *,"Brill wave + Rotating Distorted BH solve" + write(*,123)amp,eta0,sigma,byJ,mm + print*,'etamax=',etamax + 123 format(1x, 'Pars: amp',f8.5,' eta0',f8.5,' sigma',f8.5,' byJ',f9.5 + $ ,' mm',f8.5) + +c Sovle on this sized cartesian grid +c 3D grid size NE x NT +c Add 2 zones for eta coordinate and 4 for theta +c and phi coordenate. + ne = ne+2 + nq = nq+2 + np = np+2 +c + allocate(ac(ne,nq,np),ae(ne,nq,np),aw(ne,nq,np),an(ne,nq,np), + $ as(ne,nq,np),aq(ne,nq,np),ab(ne,nq,np),rhs(ne,nq,np), + $ Ksq(ne,nq,np),psisph(ne,nq,np),psip(ne,nq,np), + $ psipp(ne,nq,np),detapsisph(ne,nq,np),dqpsisph(ne,nq,np), + $ dphipsisph(ne,nq,np),detaetapsisph(ne,nq,np), + $ detaqpsisph(ne,nq,np),detaphipsisph(ne,nq,np), + $ dqqpsisph(ne,nq,np),dqphipsisph(ne,nq,np), + $ dphiphipsisph(ne,nq,np), + $ psi3d(nx,ny,nz),detapsi3d(nx,ny,nz),dqpsi3d(nx,ny,nz), + $ dphipsi3d(nx,ny,nz),detaetapsi3d(nx,ny,nz), + $ detaqpsi3d(nx,ny,nz),detaphipsi3d(nx,ny,nz), + $ dqqpsi3d(nx,ny,nz),dqphipsi3d(nx,ny,nz), + $ dphiphipsi3d(nx,ny,nz)) + allocate(etagrd(ne),qgrd(nq),phigrd(np)) + allocate(eta(nx,ny,nz),abseta(nx,ny,nz),sign_eta(nx,ny,nz), + $ q(nx,ny,nz),phi(nx,ny,nz)) + allocate(r_BL(nx,ny,nz),Delta(nx,ny,nz),Rho(nx,ny,nz)) + allocate(beta_phi(nx,ny,nz)) +c +c Initialize some array +c + nocts = 4 + nquads = 2 + dphi = nocts*0.5*pi/(np-2) + dq = nquads*0.5*pi/(nq-2) + deta = etamax/(ne-3) + + do k = 1,np + phigrd(k) = (k-1.5)*dphi + enddo + do j = 1,nq + qgrd(j) = (j-1.5)*dq + enddo + do i=1,ne + etagrd(i) = (i-2)*deta + enddo + +c +c Specify the initial guess for the conformal factor: +c + psip = 0. + psipp = 0. + do k = 1,np + do j = 1,nq + do i = 1,ne + psisph(i,j,k) = 2.*cosh(0.5*etagrd(i)) + enddo + enddo + enddo +c +c Compute K^2: +c + do k = 1,np + do j = 1,nq + do i = 1,ne +#include "CactusEinstein/RotatingDBHIVP/src/gauss.x" +#include "CactusEinstein/RotatingDBHIVP/src/ksq_odd.x" + enddo + enddo + enddo +c +c Solve the Hamiltonian constraint for the conformal factor: +c + ntries = 0 + + do it = 1,itmax+1 + + ac = 0. + ae = 0. + aw = 0. + an = 0. + as = 0. + aq = 0. + ab = 0. + rhs = 0. + +c + do k = 2,np-1 + do j = 2,nq-1 + do i = 2,ne-1 + rhs(i,j,k) = + $ (psip(i+1,j,k)-2.*psip(i,j,k)+psip(i-1,j,k))/ + $ deta**2+(psip(i,j+1,k)-2.*psip(i,j,k)+ + $ psip(i,j-1,k))/dq**2+0.5*(psip(i,j+1,k)- + $ psip(i,j-1,k))/(dq*tan(qgrd(j)))+ + $ (psip(i,j,k+1)-2.*psip(i,j,k)+psip(i,j,k-1))/ + $ (dphi**2*sin(qgrd(j))**2)-0.25*psip(i,j,k)+0.125* + $ Ksq(i,j,k)/(psisph(i,j,k)+psip(i,j,k))**7 + enddo + enddo + enddo +c + rhs = -rhs + +c +c See if the residual is small enough: +c + rhsmax = 0. +c + do k = 2,np-1 + do j = 2,nq-1 + do i = 2,ne-1 + rhsmax = max(rhsmax,abs(rhs(i,j,k))) + enddo + enddo + enddo + + if (verbose == 1) then + ntries = ntries+1 + print*,'-----------Retry', ntries + print*,'residual =', rhsmax + print*,'--------------------------------------------------' + endif + + if (rhsmax.lt.rbh_tol) goto 110 +c +c Compute stencil coefficients: +c + do k = 2,np-1 + do j = 2,nq-1 + do i = 2,ne-1 + ac(i,j,k) = -2./deta**2-2./dq**2-2./(dphi**2* + $ sin(qgrd(j))**2)-0.25-0.875*Ksq(i,j,k)/(psisph(i,j,k)+ + $ psip(i,j,k))**8 +c + ae(i,j,k) = 1./deta**2 +c + aw(i,j,k) = 1./deta**2 +c + an(i,j,k) = 1./dq**2 + 0.5/(tan(qgrd(j))*dq) +c + as(i,j,k) = 1./dq**2 - 0.5/(tan(qgrd(j))*dq) +c + aq(i,j,k) = 1./(sin(qgrd(j))**2*dphi**2) +c + ab(i,j,k) = 1./(sin(qgrd(j))**2*dphi**2) + enddo + enddo + enddo +c +c Apply boundary conditions: +c +c i=2: +c + do k = 3,np-2 + do j = 3,nq-2 + ae(2,j,k) = ae(2,j,k) + aw(2,j,k) + aw(2,j,k) = 0. +c +c i=ne-1: +c + ac(ne-1,j,k) = ac(ne-1,j,k) + 4.* + & ae(ne-1,j,k)/(3.+deta) + aw(ne-1,j,k) = aw(ne-1,j,k) - + & ae(ne-1,j,k)/(3.+deta) + ae(ne-1,j,k) = 0. + enddo +c +c j=2: +c + do i = 3,ne-2 + ac(i,2,k) = ac(i,2,k) + as(i,2,k) + as(i,2,k) = 0. +c +c j=nq-1: +c + ac(i,nq-1,k) = ac(i,nq-1,k) + + & an(i,nq-1,k) + an(i,nq-1,k) = 0. + enddo + enddo +c +c k=2: +c + do j = 3,nq-2 + do i = 3,ne-2 + ac(i,j,2) = ac(i,j,2) + ab(i,j,2) + ab(i,j,2) = 0. +c +c k=np-1: +c + ac(i,j,np-1) = ac(i,j,np-1) + + & aq(i,j,np-1) + aq(i,j,np-1) = 0. + enddo + enddo +c +c i=2, j=2: +c + do k = 3,np-2 + ae(2,2,k) = ae(2,2,k) + aw(2,2,k) + ac(2,2,k) = ac(2,2,k) + as(2,2,k) + aw(2,2,k) = 0. + as(2,2,k) = 0. +c +c i=2, j=nq-1: +c + ae(2,nq-1,k) = ae(2,nq-1,k) + aw(2,nq-1,k) + ac(2,nq-1,k) = ac(2,nq-1,k) + an(2,nq-1,k) + aw(2,nq-1,k) = 0. + an(2,nq-1,k) = 0. +c +c i=ne-1, j=2: +c + ac(ne-1,2,k) = ac(ne-1,2,k) + 4.*ae(ne-1,2,k)/(3.+deta) + + & as(ne-1,2,k) + aw(ne-1,2,k) = aw(ne-1,2,k) - ae(ne-1,2,k)/(3.+deta) + ae(ne-1,2,k) = 0. + as(ne-1,2,k) = 0. +c +c i=ne-1, j=nq-1: +c + ac(ne-1,nq-1,k) = ac(ne-1,nq-1,k) + 4.* + & ae(ne-1,nq-1,k)/(3.+deta) + an(ne-1,nq-1,k) + aw(ne-1,nq-1,k) = aw(ne-1,nq-1,k) - + & ae(ne-1,nq-1,k)/(3.+deta) + ae(ne-1,nq-1,k) = 0. + an(ne-1,nq-1,k) = 0. + enddo +c +c i=2, k=2: +c + do j = 3,nq-2 + ae(2,j,2) = ae(2,j,2) + aw(2,j,2) + ac(2,j,2) = ac(2,j,2) + ab(2,j,2) + aw(2,j,2) = 0. + ab(2,j,2) = 0. +c +c i=2, k=np-1: +c + ae(2,j,np-1) = ae(2,j,np-1) + aw(2,j,np-1) + ac(2,j,np-1) = ac(2,j,np-1) + aq(2,j,np-1) + aw(2,j,np-1) = 0. + aq(2,j,np-1) = 0. +c +c i=ne-1, k=2: +c + ac(ne-1,j,2) = ac(ne-1,j,2) + 4.*ae(ne-1,j,2)/(3.+deta) + + & ab(ne-1,j,2) + aw(ne-1,j,2) = aw(ne-1,j,2) - ae(ne-1,j,2)/(3.+deta) + ae(ne-1,j,2) = 0. + ab(ne-1,j,2) = 0. +c +c i=ne-1, k=np-1: +c + ac(ne-1,j,np-1) = ac(ne-1,j,np-1) + 4.* + & ae(ne-1,j,np-1)/(3.+deta) + aq(ne-1,j,np-1) + aw(ne-1,j,np-1) = aw(ne-1,j,np-1) - + & ae(ne-1,j,np-1)/(3.+deta) + ae(ne-1,j,np-1) = 0. + aq(ne-1,j,np-1) = 0. + enddo +c +c j=2, k=2: +c + do i = 3,ne-2 + ac(i,2,2) = ac(i,2,2) + as(i,2,2) + ab(i,2,2) + as(i,2,2) = 0. + ab(i,2,2) = 0. +c +c j=2, k=np-1: +c + ac(i,2,np-1) = ac(i,2,np-1) + as(i,2,np-1) + + & aq(i,2,np-1) + as(i,2,np-1) = 0. + aq(i,2,np-1) = 0. +c +c j=nq-1, k=2: +c + ac(i,nq-1,2) = ac(i,nq-1,2) + an(i,nq-1,2) + + & ab(i,nq-1,2) + an(i,nq-1,2) = 0. + ab(i,nq-1,2) = 0. +c +c j=nq-1, k=np-1: +c + ac(i,nq-1,np-1) = ac(i,nq-1,np-1) + an(i,nq-1,np-1) + + & aq(i,nq-1,np-1) + an(i,nq-1,np-1) = 0. + aq(i,nq-1,np-1) = 0. + enddo +c +c i=2, j=2, k=2: +c + ae(2,2,2) = ae(2,2,2) + aw(2,2,2) + ac(2,2,2) = ac(2,2,2) + as(2,2,2) + ab(2,2,2) + aw(2,2,2) = 0. + as(2,2,2) = 0. + ab(2,2,2) = 0. +c +c i=2, j=2, k=np-1: +c + ae(2,2,np-1) = ae(2,2,np-1) + aw(2,2,np-1) + ac(2,2,np-1) = ac(2,2,np-1) + as(2,2,np-1) + aq(2,2,np-1) + aw(2,2,np-1) = 0. + as(2,2,np-1) = 0. + aq(2,2,np-1) = 0. +c +c i=2, j=nq-1, k=2: +c + ae(2,nq-1,2) = ae(2,nq-1,2) + aw(2,nq-1,2) + ac(2,nq-1,2) = ac(2,nq-1,2) + an(2,nq-1,2) + ab(2,nq-1,2) + aw(2,nq-1,2) = 0. + an(2,nq-1,2) = 0. + ab(2,nq-1,2) = 0. +c +c i=2, j=nq-1, k=np-1: +c + ae(2,nq-1,np-1) = ae(2,nq-1,np-1) + aw(2,nq-1,np-1) + ac(2,nq-1,np-1) = ac(2,nq-1,np-1) + an(2,nq-1,np-1) + aq(2,nq-1,np-1) + aw(2,nq-1,np-1) = 0. + an(2,nq-1,np-1) = 0. + aq(2,nq-1,np-1) = 0. +c +c i=ne-1, j=2, k=2: +c + ac(ne-1,2,2) = ac(ne-1,2,2) + 4.*ae(ne-1,2,2)/(3.+deta) + + & as(ne-1,2,2) + ab(ne-1,2,2) + aw(ne-1,2,2) = aw(ne-1,2,2) - ae(ne-1,2,2)/(3.+deta) + ae(ne-1,2,2) = 0. + as(ne-1,2,2) = 0. + ab(ne-1,2,2) = 0. +c +c i=ne-1, j=2, k=np-1: + ac(ne-1,2,np-1) = ac(ne-1,2,np-1) + 4.*ae(ne-1,2,np-1)/(3.+deta) + + & as(ne-1,2,np-1) + aq(ne-1,2,np-1) + aw(ne-1,2,np-1) = aw(ne-1,2,np-1) - ae(ne-1,2,np-1)/(3.+deta) + ae(ne-1,2,np-1) = 0. + as(ne-1,2,np-1) = 0. + aq(ne-1,2,np-1) = 0. +c +c i=ne-1, j=nq-1, k=2: +c + ac(ne-1,nq-1,2) = ac(ne-1,nq-1,2) + 4.*ae(ne-1,nq-1,2)/(3.+deta) + + & an(ne-1,nq-1,2) + ab(ne-1,nq-1,2) + aw(ne-1,nq-1,2) = aw(ne-1,nq-1,2) - ae(ne-1,nq-1,2)/(3.+deta) + ae(ne-1,nq-1,2) = 0. + an(ne-1,nq-1,2) = 0. + ab(ne-1,nq-1,2) = 0. +c +c i=ne-1, j=nq-1, k=np-1: +c + ac(ne-1,nq-1,np-1) = ac(ne-1,nq-1,np-1) + + & 4.*ae(ne-1,nq-1,np-1)/(3.+deta) + an(ne-1,nq-1,np-1) + + & aq(ne-1,nq-1,np-1) + aw(ne-1,nq-1,np-1) = aw(ne-1,nq-1,np-1) - ae(ne-1,nq-1,np-1)/ + $ (3.+deta) + ae(ne-1,nq-1,np-1) = 0. + an(ne-1,nq-1,np-1) = 0. + aq(ne-1,nq-1,np-1) = 0. +c +c Call stab routine: +c + call rbicgst3d(ac,ae,aw,an,as,aq,ab,psipp,rhs,rbh_eps,rmax,ier, + $ ne,nq,np) + if (rmax.gt.1.0e-10) then + write(*,*) '***WARNING: bicgst3d did not converge.' + endif + if (ier.eq.-1) then + write(*,*) '***WARNING: ier=-1' + endif +c +c Now, apply boundary conditions to psipp: +c + do k = 1,np + do j = 1,nq + psipp(1,j,k) = psipp(3,j,k) + psipp(ne,j,k) = (4.*psipp(ne-1,j,k)-psipp(ne-2,j,k))/(3.+ + $ deta) + enddo + do i = 1,ne + psipp(i,1,k) = psipp(i,2,k) + psipp(i,nq,k) = psipp(i,nq-1,k) + enddo + enddo + do j = 1,nq + do i = 1,ne + psipp(i,j,1) = psipp(i,j,2) + psipp(i,j,np) = psipp(i,j,np-1) + enddo + enddo +c +c Update psip: +c + do k = 1,np + do j = 1,nq + do i = 1,ne + psip(i,j,k) = psip(i,j,k) + psipp(i,j,k) + enddo + enddo + enddo +c + enddo +c + write(*,*) 'It did not converge.' + stop +c + 110 continue + + write(*,*) '--------------------------------------------------' + print*,'Converge at Residual = ',rhsmax + +c +c Here, compute the derivatives of the spherical conformal factor +c + do k = 1,np + do j = 1,nq + do i = 2,ne-1 + eta = etagrd(i) + detapsisph(i,j,k)=0.5*(psip(i+1,j,k)-psip(i-1,j,k))/ + $ deta + sinh(0.5*etagrd(i)) + enddo + detapsisph(1,j,k) = -detapsisph(3,j,k) + enddo + enddo +c + do k = 1,np + do j = 2,nq-1 + do i = 1,ne + dqpsisph(i,j,k)=0.5*(psip(i,j+1,k)-psip(i,j-1,k))/dq + enddo + enddo + do i = 1,ne + dqpsisph(i,1,k) = -dqpsisph(i,2,k) + dqpsisph(i,nq,k) = -dqpsisph(i,nq-1,k) + enddo + enddo +c + do k = 2,np-1 + do j = 1,nq + do i = 1,ne + dphipsisph(i,j,k)=0.5*(psip(i,j,k+1)-psip(i,j,k-1))/ + $ dphi + enddo + enddo + enddo + do j = 1,nq + do i = 1,ne + dphipsisph(i,j,1) = -dphipsisph(i,j,2) + dphipsisph(i,j,np) = -dphipsisph(i,j,np-1) + enddo + enddo +c + do k = 1,np + do j = 1,nq + do i = 2,ne-1 + eta = etagrd(i) + detaetapsisph(i,j,k)=(psip(i+1,j,k)-2.*psip(i,j,k)+ + & psip(i-1,j,k))/deta**2 + sqrt(0.25)* + & cosh(0.5*etagrd(i)) + enddo + detaetapsisph(1,j,k) = detaetapsisph(3,j,k) + enddo + enddo +c + do k = 1,np + do j = 2,nq-1 + do i = 1,ne + detaqpsisph(i,j,k)=0.5*(detapsisph(i,j+1,k)- + $ detapsisph(i,j-1,k))/dq + enddo + enddo + do i = 1,ne + detaqpsisph(i,1,k) = -detaqpsisph(i,2,k) + detaqpsisph(i,nq,k) = -detaqpsisph(i,nq-1,k) + enddo + enddo +c + do k = 2,np-1 + do j = 1,nq + do i = 1,ne + detaphipsisph(i,j,k)=0.5*(detapsisph(i,j,k+1)- + $ detapsisph(i,j,k-1))/dphi + enddo + enddo + enddo + do j = 1,nq + do i = 1,ne + detaphipsisph(i,j,1) = -detaphipsisph(i,j,2) + detaphipsisph(i,j,np) = -detaphipsisph(i,j,np-1) + enddo + enddo +c + do k = 1,np + do j = 2,nq-1 + do i = 1,ne + dqqpsisph(i,j,k)=0.5*(dqpsisph(i,j+1,k)- + $ dqpsisph(i,j-1,k))/dq + enddo + enddo + do i = 1,ne + dqqpsisph(i,1,k) = dqqpsisph(i,2,k) + dqqpsisph(i,nq,k) = dqqpsisph(i,nq-1,k) + enddo + enddo +c + do k = 2,np-1 + do j = 1,nq + do i = 1,ne + dqphipsisph(i,j,k)=0.5*(dqpsisph(i,j,k+1)- + $ dqpsisph(i,j,k-1))/dphi + enddo + enddo + enddo + do j = 1,nq + do i = 1,ne + dqphipsisph(i,j,1) = -dqphipsisph(i,j,2) + dqphipsisph(i,j,np) = -dqphipsisph(i,j,np-1) + enddo + enddo +c + do k = 2,np-1 + do j = 1,nq + do i = 1,ne + dphiphipsisph(i,j,k)=0.5*(dphipsisph(i,j,k+1)- + $ dphipsisph(i,j,k-1))/dphi + enddo + enddo + enddo + do j = 1,nq + do i = 1,ne + dphiphipsisph(i,j,1) = dphiphipsisph(i,j,2) + dphiphipsisph(i,j,np) = dphiphipsisph(i,j,np-1) + enddo + enddo +c + do k = 1,np + do j = 1,nq + psisph(:,j,k)=psip(:,j,k)+2.0*cosh(0.5*etagrd) + enddo + enddo +c +c Now compute on the Cartesian coordinate. +c +c Compute eta,q,phi at the each points of cartesian grid +c + + eta = 0.5d0 * dlog (x**2 + y**2 + z**2) + abseta = abs (eta) + q = datan2 (sqrt (x**2 + y**2), z) + phi = datan2 (y, x) + + do k = 1,nz + do j = 1,ny + do i = 1,nx + + if(eta(i,j,k) .lt. 0)then + sign_eta(i,j,k) = -1.0d0 + else + sign_eta(i,j,k) = 1.0d0 + endif + + enddo + enddo + enddo + + call CCTK_GetInterpHandle (handle, "simple_local") + + npoints = nx*ny*nz + + call CCTK_Interp (ierror,cctkGH,handle,npoints,3,10,10, + $ ne,nq,np,abseta,q,phi, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ etagrd(1),qgrd(1),phigrd(1)-pi,deta,dq,dphi, + $ psisph,detapsisph,dqpsisph,dphipsisph,detaetapsisph, + $ detaqpsisph,detaphipsisph,dqqpsisph,dqphipsisph,dphiphipsisph, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ CCTK_VARIABLE_REAL, + $ psi3d,detapsi3d,dqpsi3d,dphipsi3d,detaetapsi3d,detaqpsi3d, + $ detaphipsi3d,dqqpsi3d,dqphipsi3d,dphiphipsi3d, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ CCTK_VARIABLE_REAL) + + psi = psi3d*exp(-0.5*eta) + detapsi3d = sign_eta*detapsi3d + detaqpsi3d = sign_eta*detaqpsi3d + detaphipsi3d = sign_eta*detaphipsi3d + + do k=1,nz + do j=1,ny + do i=1,nx +c psix = \partial psi / \partial x / psi +#include "CactusEinstein/RotatingDBHIVP/src/psi_1st_deriv.x" + +c psixx = \partial^2\psi / \partial x^2 / psi +#include "CactusEinstein/RotatingDBHIVP/src/psi_2nd_deriv.x" + enddo + enddo + enddo + +c Conformal metric (here, we define conformally flat) +c gxx = 1... + +c Derivatives of the metric +c dxgxx = 1/2 \partial gxx / \partial x = 0... +c + do k= 1,nz + do j= 1,ny + do i= 1,nx +#include "CactusEinstein/RotatingDBHIVP/src/gij.x" + enddo + enddo + enddo +c Extrinsic Curvture:Cactus only needs physical extrinsic curvature + do k= 1,nz + do j= 1,ny + do i= 1,nx +#include "CactusEinstein/RotatingDBHIVP/src/cgauss.x" +#include "CactusEinstein/RotatingDBHIVP/src/kij_odd.x" + enddo + enddo + enddo + + kxx = kxx/psi**2 + kxy = kxy/psi**2 + kxz = kxz/psi**2 + kyy = kyy/psi**2 + kyz = kyz/psi**2 + kzz = kzz/psi**2 + +c Set ADM mass + i = ne-15 + adm = 0.0 + do k=2,np-1 + do j=2,nq-1 + adm=adm+(psisph(i,j,k)-(psisph(i+1,j,k)-psisph(i-1,j,k))/ + $ deta)*exp(0.5*etagrd(i)) + enddo + enddo + adm=adm/(nq-2)/(np-2) + print *,'ADM mass: ',adm + if (CCTK_Equals(initial_lapse,"schwarz")==1) then + write (*,*)"Initial with schwarzschild-like lapse" + write (*,*)"using alp = (2.*r - adm)/(2.*r+adm)." + alp = (2.*r - adm)/(2.*r+adm) + endif +c + print*,'Angular momentum parameter: a/m = ', byJ/adm**2 +c +c kerr functions +c +c isotropic coordinate for schwarzschild +c + eta = 0.5d0*dlog(x**2+y**2+z**2) + $ - dlog(sqrt(adm**2-(byJ/adm)**2)/2.) + + r_BL = sqrt(adm**2-(byJ/adm)**2)/2. * exp(eta) * (1 + + $ sqrt((adm + byJ/adm)/(adm - byJ/adm)) * exp(- eta)) * (1 + + $ sqrt((adm - byJ/adm)/(adm + byJ/adm)) * exp(- eta)) + + Delta = r_BL**2 -2 * adm * r_BL + (byJ/adm)**2 + + Rho = sqrt(r_BL**2 + (byJ/adm)**2 * cos(q)**2) + +c kerr lapse for cartesian coordinate + if (kerr_slice==1) then + write (*,*) "Initial with kerr-like lapse" + alp = 1. / (1+2 * adm * r_BL * ((r_BL**2+(byJ/adm)**2)/( + $ Delta * Rho**2))) + endif + +c kerr shift for cartesian coordinate + if (kerr_shift==1) then + beta_phi = -2 * adm * r_BL/((r_BL**2+(byJ/adm)**2)**2 - + $ (byJ/adm)**2 * Delta * sin(q)**2) + + betax = -y * beta_phi + + betay = x * beta_phi + + betaz = 0.0D0 + endif + + deallocate(ac,ae,aw,an,as,aq,ab,rhs,Ksq,psisph,psip,psipp, + $ detapsisph,dqpsisph,dphipsisph,detaetapsisph,detaqpsisph, + $ detaphipsisph,dqqpsisph,dqphipsisph,dphiphipsisph, + $ psi3d,detapsi3d,dqpsi3d, + $ dphipsi3d,detaetapsi3d,detaqpsi3d, + $ detaphipsi3d,dqqpsi3d,dqphipsi3d, + $ dphiphipsi3d) + deallocate(etagrd,qgrd,phigrd) + deallocate(eta,abseta,sign_eta,q,phi) + deallocate(r_BL,Delta,Rho,beta_phi) + + return + end + + + diff --git a/src/Stab3d.F b/src/Stab3d.F new file mode 100644 index 0000000..0b1f97d --- /dev/null +++ b/src/Stab3d.F @@ -0,0 +1,277 @@ + subroutine rbicgst3d(cc,ce,cw,cn,cs,ct,cb,u,rhs, + & eps,rmax,ier,im,jm,km) +c +c This routine was lifted from stab.f. Minor modifications have +c been made. +c + implicit none +c + integer,intent(in) :: im,jm,km + real*8,intent(inout) :: cc(im,jm,km),cn(im,jm,km),cs(im,jm,km), + $ ce(im,jm,km),cw(im,jm,km),ct(im,jm,km),cb(im,jm,km) + real*8,intent(out) :: eps + real*8,intent(out) :: rmax + real*8 :: u(im,jm,km),rhs(im,jm,km) +c Local variable + integer ncyc + integer iscale,i,j,k,ier +c +c****************************************************************************** +c + + ncyc = (im-2)*(jm-2)*(km-2) + + ier=0 + +* +* Determine whether we can diagonally scale the problem to speed +* convergence. Can only be done if there are no zeros on the main +* diagonal (ie. central difference coefficient). +* + iscale=1 + do k = 2,km-1 + do j = 2,jm-1 + do i = 2,im-1 + if (cc(i,j,k).eq.0.) iscale = 0 + enddo + enddo + enddo + +* +* Do the diagonal scaling if we can. +* + if (iscale.eq.1) then + + + do k = 2,km-1 + do j = 2,jm-1 + do i = 2,im-1 + rhs(i,j,k)=rhs(i,j,k)/cc(i,j,k) + cb(i,j,k)=cb(i,j,k)/cc(i,j,k) + ct(i,j,k)=ct(i,j,k)/cc(i,j,k) + cw(i,j,k)=cw(i,j,k)/cc(i,j,k) + ce(i,j,k)=ce(i,j,k)/cc(i,j,k) + cs(i,j,k)=cs(i,j,k)/cc(i,j,k) + cn(i,j,k)=cn(i,j,k)/cc(i,j,k) + cc(i,j,k)=1. + enddo + enddo + enddo + + endif + +* +* Now call the bicgstab routine +* + + call rbicgstab (cc,cn,cs,ce,cw,ct,cb,u,rhs,eps,ncyc,rmax,ier, + & im,jm,km) + if (rmax.gt.eps) then + ier = -1 + print*,'Did not converge' + print*,' maximum residual = ',rmax + print*,' tolerance = ',eps + endif + + return + end +c +c****************************************************************** +c + subroutine rbicgstab (cc,cn,cs,ce,cw,ct,cb,x,r,tol,ncyc,rnorm, + $ ier,im,jm,km) +c +c This routine was lifted from stab.f. Minor modifications have +c been made. +c + implicit none +c + integer,intent(in) :: ncyc,im,jm,km + real*8,intent(out) :: cc(im*jm*km),cn(im*jm*km) + real*8,intent(out) :: cs(im*jm*km),ce(im*jm*km) + real*8,intent(out) :: cw(im*jm*km),ct(im*jm*km) + real*8,intent(out) :: cb(im*jm*km) + real*8,intent(out) :: tol,rnorm + integer,intent(out) :: ier + real*8 x(im*jm*km),r(im*jm*km) +c Local variables + integer :: i,j,k,kk + real*8 :: p(im*jm*km),Ap(im*jm*km),w(im*jm*km),As(im*jm*km) + real*8 :: omega, chi,chi1,chi2, delta, deltap, pp +* +*********************************************************************** +* + + do i = 1,im*jm*km + p(i) = 0. + Ap(i) = 0. + w(i) = 0. + As(i) = 0. + enddo + + kk = 0 + 10 call rusermv(cc,cn,cs,ce,cw,ct,cb,x,Ap,im,jm,km) + + do i = 1,im*jm*km + r(i) = r(i)-Ap(i) + p(i) = r(i) + enddo + +c delta = sum(r) + delta = 0. + + do i = 1,im*jm*km + delta = delta+r(i) + enddo + + if (delta.eq.0.) then + ier=-1 + return + endif + + call rusermv(cc,cn,cs,ce,cw,ct,cb,p,Ap,im,jm,km) + +c phi = sum(Ap) + pp = 0. + + do i = 1,im*jm*km + pp = pp+Ap(i) + enddo + pp = pp/delta + + if (pp.eq.0.) then + ier=-1 + return + endif + +c rnorm = sum(r**2) + rnorm = 0. + + do i = 1,im*jm*km + rnorm = rnorm + r(i)**2 + enddo + rnorm=sqrt(rnorm) + +c Test if initial guess is great (residual less than tolerance) + if (rnorm .lt. tol) return + + 1 continue + kk = kk + 1 + + omega = 1./pp + + + do i = 1,im*jm*km + w(i) = r(i) - omega*Ap(i) + enddo + + call rusermv(cc,cn,cs,ce,cw,ct,cb,w,As,im,jm,km) + +c chi1 = sum(As*w) + chi1 = 0. + + do i = 1,im*jm*km + chi1 = chi1+As(i)*w(i) + enddo +c chi2 = sum(As**2) + chi2 = 0. + + do i = 1,im*jm*km + chi2 = chi2+As(i)**2 + enddo + + chi = chi1/chi2 + + + do i = 1,im*jm*km + r(i) = w(i) - chi*As(i) + x(i) = x(i) + omega*p(i) + chi*w(i) + enddo + + deltap = delta +c delta = sum(r) + delta = 0. + + do i = 1,im*jm*km + delta = delta+r(i) + enddo + + if (delta.eq.0.) then + goto 10 + endif + + + do i = 1,im*jm*km + p(i) = r(i) + (p(i)-chi*Ap(i))*omega* + & delta/(deltap*chi) + enddo + + call rusermv(cc,cn,cs,ce,cw,ct,cb,p,Ap,im,jm,km) + +c phi = sum(Ap) + pp = 0. + + do i = 1,im*jm*km + pp = pp+Ap(i) + enddo + pp=pp/delta + + if (pp.eq.0.) then + goto 10 + endif + + if (kk.gt.ncyc) then + print*,' BI-CGStab solver reached maximum nuber of iterations.' + ier=-1 + return + endif + +c rnorm = sum(r**2) + rnorm = 0. + + do i = 1,im*jm*km + rnorm = rnorm+r(i)**2 + enddo + rnorm=sqrt(rnorm) + + if (rnorm .gt. tol) goto 1 + write(*,*) "ncyc=",kk + return + end +c +c****************************************************************** +c + subroutine rusermv(cc,cn,cs,ce,cw,ct,cb,x,y,im,jm,km) +c +c This routine was lifted from stab.f. Minor modifications have +c been made. +c +c Be careful that the c's are zero on their outer boundary!! +c + implicit none +c + integer,intent(in) :: im,jm,km + real*8,intent(out) :: cc(im*jm*km),cn(im*jm*km) + real*8,intent(out) :: cs(im*jm*km),ce(im*jm*km) + real*8,intent(out) :: cw(im*jm*km),ct(im*jm*km) + real*8,intent(out) :: cb(im*jm*km) + real*8 :: x(im*jm*km), y(im*jm*km) +c Local variables + integer :: i, j, k +* +*********************************************************************** +* +* + do i = im*jm+im+1,im*(jm*km-jm-1) + y(i) = cw(i)*x(i-1) + & +cc(i)*x(i) + & +ce(i)*x(i+1) + & +cn(i)*x(i+im) + & +cs(i)*x(i-im) + & +ct(i)*x(i+im*jm) + & +cb(i)*x(i-im*jm) + enddo +c + return + end + diff --git a/src/cgauss.x b/src/cgauss.x new file mode 100644 index 0000000..dda383d --- /dev/null +++ b/src/cgauss.x @@ -0,0 +1,57 @@ + gtil = amp*(exp(-((-eta0 + eta(i,j,k))**2/sigma**2)) + exp(-((eta0 + + & eta(i,j,k))**2/sigma**2))) + dngtil = amp*((-2.00000000000000d0*(-eta0 + eta(i,j,k))*exp(-((-eta0 + & + eta(i,j,k))**2/sigma**2)))/sigma**2 - (2.00000000000000d0*(eta0 + + & eta(i,j,k))*exp(-((eta0 + eta(i,j,k))**2/sigma**2)))/sigma**2) + dnngtil = amp*((4.0000000000000d0*(-eta0 + eta(i,j,k))**2*exp(-((-et + & a0 + eta(i,j,k))**2/sigma**2)))/sigma**4 - (2.00000000000000d0*exp(- + & ((-eta0 + eta(i,j,k))**2/sigma**2)))/sigma**2 + (4.0000000000000d0*( + & eta0 + eta(i,j,k))**2*exp(-((eta0 + eta(i,j,k))**2/sigma**2)))/sigma**4 + & - (2.00000000000000d0*exp(-((eta0 + eta(i,j,k))**2/sigma**2)))/sigma + & **2) + dnnngtil = amp*((-8.0000000000000d0*(-eta0 + eta(i,j,k))**3*exp(-((- + & eta0 + eta(i,j,k))**2/sigma**2)))/sigma**6 + (1.20000000000000d1*(-e + & ta0 + eta(i,j,k))*exp(-((-eta0 + eta(i,j,k))**2/sigma**2)))/sigma**4 - ( + & 8.0000000000000d0*(eta0 + eta(i,j,k))**3*exp(-((eta0 + eta(i,j,k))**2/si + & gma**2)))/sigma**6 + (1.20000000000000d1*(eta0 + eta(i,j,k))*exp(-(( + & eta0 + eta(i,j,k))**2/sigma**2)))/sigma**4) + dnnnngtil = amp*((1.60000000000000d1*(-eta0 + eta(i,j,k))**4*exp(-(( + & -eta0 + eta(i,j,k))**2/sigma**2)))/sigma**8 - (4.8000000000000d1*(-e + & ta0 + eta(i,j,k))**2*exp(-((-eta0 + eta(i,j,k))**2/sigma**2)))/sigma**6 + & + (1.20000000000000d1*exp(-((-eta0 + eta(i,j,k))**2/sigma**2)))/sigm + & a**4 + (1.60000000000000d1*(eta0 + eta(i,j,k))**4*exp(-((eta0 + eta(i,j,k) + & )**2/sigma**2)))/sigma**8 - (4.8000000000000d1*(eta0 + eta(i,j,k)) + & **2*exp(-((eta0 + eta(i,j,k))**2/sigma**2)))/sigma**6 + (1.200000000 + & 00000d1*exp(-((eta0 + eta(i,j,k))**2/sigma**2)))/sigma**4) + dnnnnngtil = amp*((-3.2000000000000d1*(-eta0 + eta(i,j,k))**5*exp(-( + & (-eta0 + eta(i,j,k))**2/sigma**2)))/sigma**10 + (1.60000000000000d2* + & (-eta0 + eta(i,j,k))**3*exp(-((-eta0 + eta(i,j,k))**2/sigma**2)))/sigma* + & *8 - (1.20000000000000d2*(-eta0 + eta(i,j,k))*exp(-((-eta0 + eta(i,j,k)) + & **2/sigma**2)))/sigma**6 - (3.2000000000000d1*(eta0 + eta(i,j,k))**5 + & *exp(-((eta0 + eta(i,j,k))**2/sigma**2)))/sigma**10 + (1.60000000000 + & 000d2*(eta0 + eta(i,j,k))**3*exp(-((eta0 + eta(i,j,k))**2/sigma**2)))/si + & gma**8 - (1.20000000000000d2*(eta0 + eta(i,j,k))*exp(-((eta0 + eta(i,j,k) + & )**2/sigma**2)))/sigma**6) + dnnnnnngtil = amp*((6.4000000000000d1*(-eta0 + eta(i,j,k))**6*exp(-( + & (-eta0 + eta(i,j,k))**2/sigma**2)))/sigma**12 - (4.8000000000000d2*( + & -eta0 + eta(i,j,k))**4*exp(-((-eta0 + eta(i,j,k))**2/sigma**2)))/sigma** + & 10 + (7.2000000000000d2*(-eta0 + eta(i,j,k))**2*exp(-((-eta0 + eta(i,j,k) + & )**2/sigma**2)))/sigma**8 - (1.20000000000000d2*exp(-((-eta0 + + & eta(i,j,k))**2/sigma**2)))/sigma**6 + (6.4000000000000d1*(eta0 + eta(i,j,k) + & )**6*exp(-((eta0 + eta(i,j,k))**2/sigma**2)))/sigma**12 - (4.8000 + & 000000000d2*(eta0 + eta(i,j,k))**4*exp(-((eta0 + eta(i,j,k))**2/sigma**2 + & )))/sigma**10 + (7.2000000000000d2*(eta0 + eta(i,j,k))**2*exp(-((eta + & 0 + eta(i,j,k))**2/sigma**2)))/sigma**8 - (1.20000000000000d2*exp(-( + & (eta0 + eta(i,j,k))**2/sigma**2)))/sigma**6) + dnnnnnnngtil = amp*((-1.28000000000000d2*(-eta0 + eta(i,j,k))**7*exp + & (-((-eta0 + eta(i,j,k))**2/sigma**2)))/sigma**14 + (1.34400000000000 + & d3*(-eta0 + eta(i,j,k))**5*exp(-((-eta0 + eta(i,j,k))**2/sigma**2)))/sig + & ma**12 - (3.3600000000000d3*(-eta0 + eta(i,j,k))**3*exp(-((-eta0 + e + & ta(i,j,k))**2/sigma**2)))/sigma**10 + (1.68000000000000d3*(-eta0 + e + & ta(i,j,k))*exp(-((-eta0 + eta(i,j,k))**2/sigma**2)))/sigma**8 - (1.28000 + & 000000000d2*(eta0 + eta(i,j,k))**7*exp(-((eta0 + eta(i,j,k))**2/sigma**2 + & )))/sigma**14 + (1.34400000000000d3*(eta0 + eta(i,j,k))**5*exp(-((et + & a0 + eta(i,j,k))**2/sigma**2)))/sigma**12 - (3.3600000000000d3*(eta0 + & + eta(i,j,k))**3*exp(-((eta0 + eta(i,j,k))**2/sigma**2)))/sigma**10 + ( + & 1.68000000000000d3*(eta0 + eta(i,j,k))*exp(-((eta0 + eta(i,j,k))**2/sigm + & a**2)))/sigma**8) diff --git a/src/gauss.x b/src/gauss.x new file mode 100644 index 0000000..a853201 --- /dev/null +++ b/src/gauss.x @@ -0,0 +1,57 @@ + gtil = amp*(exp(-((-eta0 + etagrd(i))**2/sigma**2)) + exp(-((eta0 + + & etagrd(i))**2/sigma**2))) + dngtil = amp*((-2.00000000000000d0*(-eta0 + etagrd(i))*exp(-((-eta0 + & + etagrd(i))**2/sigma**2)))/sigma**2 - (2.00000000000000d0*(eta0 + + & etagrd(i))*exp(-((eta0 + etagrd(i))**2/sigma**2)))/sigma**2) + dnngtil = amp*((4.0000000000000d0*(-eta0 + etagrd(i))**2*exp(-((-et + & a0 + etagrd(i))**2/sigma**2)))/sigma**4 - (2.00000000000000d0*exp(- + & ((-eta0 + etagrd(i))**2/sigma**2)))/sigma**2 + (4.0000000000000d0*( + & eta0 + etagrd(i))**2*exp(-((eta0 + etagrd(i))**2/sigma**2)))/sigma**4 + & - (2.00000000000000d0*exp(-((eta0 + etagrd(i))**2/sigma**2)))/sigma + & **2) + dnnngtil = amp*((-8.0000000000000d0*(-eta0 + etagrd(i))**3*exp(-((- + & eta0 + etagrd(i))**2/sigma**2)))/sigma**6 + (1.20000000000000d1*(-e + & ta0 + etagrd(i))*exp(-((-eta0 + etagrd(i))**2/sigma**2)))/sigma**4 - ( + & 8.0000000000000d0*(eta0 + etagrd(i))**3*exp(-((eta0 + etagrd(i))**2/si + & gma**2)))/sigma**6 + (1.20000000000000d1*(eta0 + etagrd(i))*exp(-(( + & eta0 + etagrd(i))**2/sigma**2)))/sigma**4) + dnnnngtil = amp*((1.60000000000000d1*(-eta0 + etagrd(i))**4*exp(-(( + & -eta0 + etagrd(i))**2/sigma**2)))/sigma**8 - (4.8000000000000d1*(-e + & ta0 + etagrd(i))**2*exp(-((-eta0 + etagrd(i))**2/sigma**2)))/sigma**6 + & + (1.20000000000000d1*exp(-((-eta0 + etagrd(i))**2/sigma**2)))/sigm + & a**4 + (1.60000000000000d1*(eta0 + etagrd(i))**4*exp(-((eta0 + etag + & rd(i))**2/sigma**2)))/sigma**8 - (4.8000000000000d1*(eta0 + etagrd(i)) + & **2*exp(-((eta0 + etagrd(i))**2/sigma**2)))/sigma**6 + (1.200000000 + & 00000d1*exp(-((eta0 + etagrd(i))**2/sigma**2)))/sigma**4) + dnnnnngtil = amp*((-3.2000000000000d1*(-eta0 + etagrd(i))**5*exp(-( + & (-eta0 + etagrd(i))**2/sigma**2)))/sigma**10 + (1.60000000000000d2* + & (-eta0 + etagrd(i))**3*exp(-((-eta0 + etagrd(i))**2/sigma**2)))/sigma* + & *8 - (1.20000000000000d2*(-eta0 + etagrd(i))*exp(-((-eta0 + etagrd(i)) + & **2/sigma**2)))/sigma**6 - (3.2000000000000d1*(eta0 + etagrd(i))**5 + & *exp(-((eta0 + etagrd(i))**2/sigma**2)))/sigma**10 + (1.60000000000 + & 000d2*(eta0 + etagrd(i))**3*exp(-((eta0 + etagrd(i))**2/sigma**2)))/si + & gma**8 - (1.20000000000000d2*(eta0 + etagrd(i))*exp(-((eta0 + etagr + & d(i))**2/sigma**2)))/sigma**6) + dnnnnnngtil = amp*((6.4000000000000d1*(-eta0 + etagrd(i))**6*exp(-( + & (-eta0 + etagrd(i))**2/sigma**2)))/sigma**12 - (4.8000000000000d2*( + & -eta0 + etagrd(i))**4*exp(-((-eta0 + etagrd(i))**2/sigma**2)))/sigma** + & 10 + (7.2000000000000d2*(-eta0 + etagrd(i))**2*exp(-((-eta0 + etagr + & d(i))**2/sigma**2)))/sigma**8 - (1.20000000000000d2*exp(-((-eta0 + + & etagrd(i))**2/sigma**2)))/sigma**6 + (6.4000000000000d1*(eta0 + eta + & grd(i))**6*exp(-((eta0 + etagrd(i))**2/sigma**2)))/sigma**12 - (4.8000 + & 000000000d2*(eta0 + etagrd(i))**4*exp(-((eta0 + etagrd(i))**2/sigma**2 + & )))/sigma**10 + (7.2000000000000d2*(eta0 + etagrd(i))**2*exp(-((eta + & 0 + etagrd(i))**2/sigma**2)))/sigma**8 - (1.20000000000000d2*exp(-( + & (eta0 + etagrd(i))**2/sigma**2)))/sigma**6) + dnnnnnnngtil = amp*((-1.28000000000000d2*(-eta0 + etagrd(i))**7*exp + & (-((-eta0 + etagrd(i))**2/sigma**2)))/sigma**14 + (1.34400000000000 + & d3*(-eta0 + etagrd(i))**5*exp(-((-eta0 + etagrd(i))**2/sigma**2)))/sig + & ma**12 - (3.3600000000000d3*(-eta0 + etagrd(i))**3*exp(-((-eta0 + e + & tagrd(i))**2/sigma**2)))/sigma**10 + (1.68000000000000d3*(-eta0 + e + & tagrd(i))*exp(-((-eta0 + etagrd(i))**2/sigma**2)))/sigma**8 - (1.28000 + & 000000000d2*(eta0 + etagrd(i))**7*exp(-((eta0 + etagrd(i))**2/sigma**2 + & )))/sigma**14 + (1.34400000000000d3*(eta0 + etagrd(i))**5*exp(-((et + & a0 + etagrd(i))**2/sigma**2)))/sigma**12 - (3.3600000000000d3*(eta0 + & + etagrd(i))**3*exp(-((eta0 + etagrd(i))**2/sigma**2)))/sigma**10 + ( + & 1.68000000000000d3*(eta0 + etagrd(i))*exp(-((eta0 + etagrd(i))**2/sigm + & a**2)))/sigma**8) diff --git a/src/gij.x b/src/gij.x new file mode 100644 index 0000000..693b57d --- /dev/null +++ b/src/gij.x @@ -0,0 +1,24 @@ + gxx(i,j,k) = 1.00000000000000d0 +c dxgxx(i,j,k) = 0 +c dygxx(i,j,k) = 0 +c dzgxx(i,j,k) = 0 + gxy(i,j,k) = 0 +c dxgxy(i,j,k) = 0 +c dygxy(i,j,k) = 0 +c dzgxy(i,j,k) = 0 + gxz(i,j,k) = 0 +c dxgxz(i,j,k) = 0 +c dygxz(i,j,k) = 0 +c dzgxz(i,j,k) = 0 + gyy(i,j,k) = 1.00000000000000d0 +c dxgyy(i,j,k) = 0 +c dygyy(i,j,k) = 0 +c dzgyy(i,j,k) = 0 + gyz(i,j,k) = 0 +c dxgyz(i,j,k) = 0 +c dygyz(i,j,k) = 0 +c dzgyz(i,j,k) = 0 + gzz(i,j,k) = 1.00000000000000d0 +c dxgzz(i,j,k) = 0 +c dygzz(i,j,k) = 0 +c dzgzz(i,j,k) = 0 diff --git a/src/kij_axi.x b/src/kij_axi.x new file mode 100644 index 0000000..e91a8e0 --- /dev/null +++ b/src/kij_axi.x @@ -0,0 +1,62 @@ + o1 = cos(phi) + o2 = cos(q) + o3 = o2**2 + o4 = -3.00000000000000d0*eta + o5 = exp(o4) + o6 = sin(phi) + o7 = sin(q) + o8 = o7**2 + o9 = dnnnnngtil + o10 = -3.6000000000000d1*o9 + o11 = dnnnnnnngtil + o12 = -4.0000000000000d0*o11*o8 + o13 = 1.16000000000000d2*o8*o9 + o14 = dnnngtil + o15 = -4.1500000000000d2*o8 + o16 = 1.50000000000000d2 + o15 + o17 = o14*o16 + o18 = dngtil + o19 = 5.0000000000000d0*o8 + o20 = -2.00000000000000d0 + o19 + o21 = 7.5000000000000d1*o18*o20 + o22 = o10 + o12 + o13 + o17 + o21 + o23 = dnnnnnngtil + o24 = 3.00000000000000d0*o23 + o25 = 2.00000000000000d0*q + o26 = cos(o25) + o27 = 4.0000000000000d0*o23*o26 + o28 = 4.0000000000000d0*q + o29 = cos(o28) + o30 = -7.0000000000000d0*o23*o29 + o31 = dnnnngtil + o32 = 2.10000000000000d1*o31 + o33 = 6.4000000000000d1*o26*o31 + o34 = 2.03000000000000d2*o29*o31 + o35 = o24 + o27 + o30 + o32 + o33 + o34 + o36 = -4.0000000000000d0*o35 + o37 = 2.68000000000000d2*o26 + o38 = 5.8100000000000d2*o29 + o39 = 1.11000000000000d2 + o37 + o38 + o40 = dnngtil + o41 = 5.0000000000000d0*o39*o40 + o42 = 2.00000000000000d1*o26 + o43 = 3.5000000000000d1*o29 + o44 = 9.0000000000000d0 + o42 + o43 + o45 = gtil + o46 = -7.5000000000000d1*o44*o45 + o47 = o36 + o41 + o46 + o48 = o1**2 + o49 = o6**2 + o50 = o7*o8 + hxx(i,j,k) = 6.6666666666667d-1*o1*o22*o3*o5*o6*o8 - 8.333333333 + & 3333d-2*o1*o47*o5*o6*o8 + hxy(i,j,k) = -3.3333333333333d-1*o22*o3*o48*o5*o8 + 4.1666666666 + & 667d-2*o47*o48*o5*o8 + 3.3333333333333d-1*o22*o3*o49*o5*o8 - 4.1 + & 666666666667d-2*o47*o49*o5*o8 + hxz(i,j,k) = -3.3333333333333d-1*o2*o22*o5*o50*o6 - 4.1666666666 + & 667d-2*o2*o47*o5*o6*o7 + hyy(i,j,k) = -6.6666666666667d-1*o1*o22*o3*o5*o6*o8 + 8.33333333 + & 33333d-2*o1*o47*o5*o6*o8 + hyz(i,j,k) = 3.3333333333333d-1*o1*o2*o22*o5*o50 + 4.16666666666 + & 67d-2*o1*o2*o47*o5*o7 + hzz(i,j,k) = 0 diff --git a/src/kij_odd.x b/src/kij_odd.x new file mode 100644 index 0000000..f9036f2 --- /dev/null +++ b/src/kij_odd.x @@ -0,0 +1,153 @@ + o1 = mm**2 + o2 = 5.0000000000000d1*o1 + o3 = -8.0000000000000d1 + o2 + o4 = 1/o3 + o5 = 2.00000000000000d0*q(i,j,k) + o6 = cos(o5) + o7 = 5.0000000000000d0*o6 + o8 = 3.00000000000000d0 + o7 + o9 = dnnnnngtil + o10 = dnnngtil + o11 = -2.60000000000000d1*o10 + o12 = dngtil + o13 = 2.50000000000000d1*o12 + o14 = o11 + o13 + o9 + o15 = -3.00000000000000d0*eta(i,j,k) + o16 = exp(o15) + o17 = sin(phi(i,j,k)) + o18 = o17**2 + o19 = mm*phi(i,j,k) + o20 = sin(o19) + o21 = sin(q(i,j,k)) + o22 = o21**2 + o23 = 2.50000000000000d1*o1 + o24 = -4.0000000000000d1 + o23 + o25 = 1/o24 + o26 = cos(phi(i,j,k)) + o27 = o26**2 + o28 = -1.60000000000000d1 + o1 + o29 = 3.00000000000000d0*o6 + o30 = 2.00000000000000d0 + o29 + o31 = 2.00000000000000d0*o28*o30*o9 + o32 = -1.07000000000000d2*o1 + o33 = 6.0900000000000d2*o6 + o34 = -1.59000000000000d2*o1*o6 + o35 = 4.0700000000000d2 + o32 + o33 + o34 + o36 = o10*o35 + o37 = 3.5000000000000d1*o1 + o38 = 5.0000000000000d0*o1 + o39 = -1.30000000000000d1 + o38 + o40 = 9.0000000000000d0*o39*o6 + o41 = -8.3000000000000d1 + o37 + o40 + o42 = 5.0000000000000d0*o12*o41 + o43 = o31 + o36 + o42 + o44 = o22**2 + o45 = cos(q(i,j,k)) + o46 = o45**2 + o47 = -3.00000000000000d0*mm*o14*o20*o22*o4*o8 + o48 = -(mm*o20*o22*o25*o43) + o49 = o47 + o48 + o50 = -8.0000000000000d0 + o38 + o51 = 1/o50 + o52 = cos(o19) + o53 = -9.0000000000000d0*o1 + o54 = 5.0000000000000d0*o50*o6 + o55 = 8.0000000000000d0 + o53 + o54 + o56 = -7.5000000000000d1*o12*o55 + o57 = 2.88000000000000d2*o9 + o58 = -1.20000000000000d1*o1*o9 + o59 = -9.4000000000000d1*o1 + o60 = -4.1500000000000d2*o6 + o61 = 1.75000000000000d2*o1*o6 + o62 = 1.15000000000000d2 + o59 + o60 + o61 + o63 = 4.0000000000000d0*o10*o62 + o64 = dnnnnnnngtil + o65 = 3.2000000000000d1*o22*o64 + o66 = -2.00000000000000d0*o1*o22*o64 + o67 = -9.2800000000000d2*o22*o9 + o68 = 1.03000000000000d2*o1*o22*o9 + o69 = o57 + o58 + o63 + o65 + o66 + o67 + o68 + o70 = 2.00000000000000d0*o69 + o71 = o56 + o70 + o72 = sin(o5) + o73 = dnnnnnngtil + o74 = 3.2000000000000d1*o73 + o75 = -2.00000000000000d0*o1*o73 + o76 = -3.2000000000000d1*o6*o73 + o77 = 2.00000000000000d0*o1*o6*o73 + o78 = dnnnngtil + o79 = -3.4000000000000d2*o78 + o80 = 5.5000000000000d1*o1*o78 + o81 = -8.1200000000000d2*o6*o78 + o82 = 1.70000000000000d1*o1*o6*o78 + o83 = 7.0000000000000d0*o1 + o84 = -2.09000000000000d2*o6 + o85 = 8.9000000000000d1*o1*o6 + o86 = -3.10000000000000d1 + o83 + o84 + o85 + o87 = dnngtil + o88 = -2.00000000000000d1*o86*o87 + o89 = 7.0000000000000d0*o6 + o90 = 1.00000000000000d0 + o89 + o91 = gtil + o92 = 7.5000000000000d1*o50*o90*o91 + o93 = o74 + o75 + o76 + o77 + o79 + o80 + o81 + o82 + o88 + o92 + o94 = -2.40000000000000d2*byJ + o95 = 1.50000000000000d2*byJ*o1 + o96 = -1.15200000000000d3*o46*o52*o78 + o97 = 7.2000000000000d1*o1*o46*o52*o78 + o98 = 3.9000000000000d1*o1 + o99 = -2.68000000000000d2*o6 + o100 = 1.00000000000000d2*o1*o6 + o101 = 4.0000000000000d0*q(i,j,k) + o102 = cos(o101) + o103 = -5.8100000000000d2*o102 + o104 = 2.45000000000000d2*o1*o102 + o105 = -1.11000000000000d2 + o100 + o103 + o104 + o98 + o99 + o106 = -5.0000000000000d0*o105*o52*o87 + o107 = -1.92000000000000d2*o22*o46*o52*o73 + o108 = 1.20000000000000d1*o1*o22*o46*o52*o73 + o109 = 2.88000000000000d2*o22*o52*o78 + o110 = -1.80000000000000d1*o1*o22*o52*o78 + o111 = 5.5680000000000d3*o22*o46*o52*o78 + o112 = -6.1800000000000d2*o1*o22*o46*o52*o78 + o113 = 3.2000000000000d1*o44*o52*o73 + o114 = -2.00000000000000d0*o1*o44*o52*o73 + o115 = -9.2800000000000d2*o44*o52*o78 + o116 = 1.03000000000000d2*o1*o44*o52*o78 + o117 = o106 + o107 + o108 + o109 + o110 + o111 + o112 + o113 + o + & 114 + o115 + o116 + o94 + o95 + o96 + o97 + o118 = 2.00000000000000d0*o117 + o119 = 2.00000000000000d1*o6 + o120 = 3.5000000000000d1*o102 + o121 = 9.0000000000000d0 + o119 + o120 + o122 = 1.87500000000000d1*o121*o50*o52*o91 + o123 = o118 + o122 + o124 = o21*o22 + kxx(i,j,k) = mm*o16*o20*o25*o27*o43*o44 + o16*o27*o46*o49 - 1.00 + & 000000000000d-1*o123*o16*o17*o22*o26*o51 - 1.00000000000000d-1*o + & 16*o17*o22*o26*o46*o51*o52*o71 + 3.00000000000000d0*mm*o14*o16*o + & 18*o20*o22*o4*o8 + 5.0000000000000d-2*mm*o16*o20*o21*o27*o45*o51 + & *o72*o93 + kxy(i,j,k) = mm*o16*o17*o20*o25*o26*o43*o44 + o16*o17*o26*o46*o4 + & 9 - 5.0000000000000d-2*o123*o16*o18*o22*o51 + 5.0000000000000d-2 + & *o123*o16*o22*o27*o51 - 5.0000000000000d-2*o16*o18*o22*o46*o51*o + & 52*o71 + 5.0000000000000d-2*o16*o22*o27*o46*o51*o52*o71 - 3.0000 + & 0000000000d0*mm*o14*o16*o17*o20*o22*o26*o4*o8 + 5.0000000000000d + & -2*mm*o16*o17*o20*o21*o26*o45*o51*o72*o93 + kxz(i,j,k) = mm*o124*o16*o20*o25*o26*o43*o45 - o16*o21*o26*o45*o + & 49 - 5.0000000000000d-2*o123*o16*o17*o21*o45*o51 + 5.00000000000 + & 00d-2*o124*o16*o17*o45*o51*o52*o71 - 2.50000000000000d-2*mm*o16* + & o20*o22*o26*o51*o72*o93 + 2.50000000000000d-2*mm*o16*o20*o26*o46 + & *o51*o72*o93 + kyy(i,j,k) = mm*o16*o18*o20*o25*o43*o44 + o16*o18*o46*o49 + 1.00 + & 000000000000d-1*o123*o16*o17*o22*o26*o51 + 1.00000000000000d-1*o + & 16*o17*o22*o26*o46*o51*o52*o71 + 3.00000000000000d0*mm*o14*o16*o + & 20*o22*o27*o4*o8 + 5.0000000000000d-2*mm*o16*o18*o20*o21*o45*o51 + & *o72*o93 + kyz(i,j,k) = mm*o124*o16*o17*o20*o25*o43*o45 - o16*o17*o21*o45*o + & 49 + 5.0000000000000d-2*o123*o16*o21*o26*o45*o51 - 5.00000000000 + & 00d-2*o124*o16*o26*o45*o51*o52*o71 - 2.50000000000000d-2*mm*o16* + & o17*o20*o22*o51*o72*o93 + 2.50000000000000d-2*mm*o16*o17*o20*o46 + & *o51*o72*o93 + kzz(i,j,k) = mm*o16*o20*o22*o25*o43*o46 + o16*o22*o49 - 5.000000 + & 0000000d-2*mm*o16*o20*o21*o45*o51*o72*o93 diff --git a/src/ksq2d_odd.x b/src/ksq2d_odd.x new file mode 100644 index 0000000..39a175f --- /dev/null +++ b/src/ksq2d_odd.x @@ -0,0 +1,23 @@ + o1 = cos(q) + o2 = 1. + o3 = o2**2 + o4 = 1/o3 + o5 = sin(q) + o6 = 2.00000000000000d0 + 3.0d0 + o7 = o5**2 + o8 = -(o6*o7) + o9 = 1.00000000000000d0 + 3.0d0 + o8 + o10 = o3**2 + o11 = o10**2 + o12 = dngtil**2 + o13 = o1**2 + o14 = gtil**2 + o15 = o9**2 +c exc33(i,j,k) = 0 +c exc32(i,j,k) = -(dngtil*o1*o4*o5**n) +c exc31(i,j,k) = gtil*o4*o5**(-1.00000000000000d0 + n)*o9 +c exc22(i,j,k) = 0 +c exc21(i,j,k) = 0 +c exc11(i,j,k) = 0 + ksq(i,j,k) = (2.00000000000000d0*o5**(2.00000000000000d0*(-2.000 + & 00000000000d0 + 3.d0))*(o14*o15 + o12*o13*o7))/(o10*o11) diff --git a/src/ksq_axi.x b/src/ksq_axi.x new file mode 100644 index 0000000..57ad7ec --- /dev/null +++ b/src/ksq_axi.x @@ -0,0 +1,111 @@ + o1 = cos(q) + o2 = 2.00000000000000d0*q + o3 = cos(o2) + o4 = dnnnnngtil + o5 = sin(q) + o6 = o5**2 + o7 = o5*o6 + o8 = dnnngtil + o9 = dngtil + o10 = dnnnnnnngtil + o11 = o6**2 + o12 = 4.0000000000000d0*q + o13 = cos(o12) + o14 = dnnnngtil + o15 = dnngtil + o16 = o1**2 + o17 = dnnnnnngtil + o18 = o11*o6 + o19 = gtil + o20 = byJ**2 + o21 = o14**2 + o22 = o3**2 + o23 = o13**2 + o24 = o15**2 + o25 = o4**2 + o26 = o8**2 + o27 = o9**2 + o28 = o16**2 + o29 = o17**2 + o30 = o10**2 + o31 = o11**2 + o32 = o31*o6 + o33 = o19**2 +c exc33(i,j,k) = 0 +c exc32(i,j,k) = 1.33333333333333d0*o1*o10*o11*o5 - 2.083333333333 +c & 33d-2*o1*(3.5200000000000d2 - 9.2800000000000d2*o3)*o4*o7 + 1.66 +c & 666666666667d-1*o1*(1.15000000000000d2 - 4.1500000000000d2*o3)*o +c & 7*o8 - 1.56250000000000d0*o1*(8.0000000000000d0 - 4.000000000000 +c & 0d1*o3)*o7*o9 +c exc31(i,j,k) = -8.0000000000000d0*o11*o16*o17 + 1.33333333333333 +c & d0*o17*o18 + 3.00000000000000d0*byJ*o6 - 2.08333333333333d-1*o15*( +c & -1.11000000000000d2 - 5.8100000000000d2*o13 - 2.68000000000000d2 +c & *o3)*o6 - 5.2083333333333d-3*o14*(6.7200000000000d2 + 6.49600000 +c & 00000d3*o13 + 2.04800000000000d3*o3)*o6 + o19*(5.0000000000000d1 +c & *o11 + 7.5000000000000d2*o11*o16 - 1.25000000000000d2*o18 - 2.00 +c & 000000000000d2*o16*o6) +c exc22(i,j,k) = 0 +c exc21(i,j,k) = 0 +c exc11(i,j,k) = 0 + t1 = -9.6000000000000d1*byJ*o11*o16*o17 + 1.12000000000000d2*o11*o + & 14*o16*o17 + 1.08266666666667d3*o11*o13*o14*o16*o17 - 7.40000000 + & 00000d2*o11*o15*o16*o17 - 3.8733333333333d3*o11*o13*o15*o16*o17 + & + 1.60000000000000d1*byJ*o17*o18 - 1.86666666666667d1*o14*o17*o18 + & - 1.80444444444444d2*o13*o14*o17*o18 + 1.23333333333333d2*o15*o1 + & 7*o18 + 6.4555555555556d2*o13*o15*o17*o18 + 6.0000000000000d2*byJ* + & o11*o19 - 7.0000000000000d2*o11*o14*o19 - 6.7666666666667d3*o11* + & o13*o14*o19 + 4.6250000000000d3*o11*o15*o19 + 2.42083333333333d4 + & *o11*o13*o15*o19 + 9.0000000000000d3*byJ*o11*o16*o19 - 1.050000000 + & 00000d4*o11*o14*o16*o19 - 1.01500000000000d5*o11*o13*o14*o16*o19 + & + 6.9375000000000d4*o11*o15*o16*o19 + 3.6312500000000d5*o11*o13 + & *o15*o16*o19 - 1.50000000000000d3*byJ*o18*o19 + 1.75000000000000d3 + & *o14*o18*o19 + 1.69166666666667d4*o13*o14*o18*o19 - 1.1562500000 + & 0000d4*o15*o18*o19 - 6.0520833333333d4*o13*o15*o18*o19 - 2.66666 + & 666666667d3*o16*o17*o18*o19 + 1.07555555555556d2*o11*o16*o25 + t2 = 7.4755555555556d2*o11*o16*o22*o25 + 7.3472222222222d2*o11*o + & 16*o26 + 9.5680555555556d3*o11*o16*o22*o26 + 3.12500000000000d2* + & o11*o16*o27 + 7.8125000000000d3*o11*o16*o22*o27 + 6.400000000000 + & 0d3*o11*o17*o19*o28 - 2.40000000000000d4*o17*o18*o19*o28 + 1.280 + & 00000000000d2*o18*o28*o29 + 3.4133333333333d2*o11*o14*o16*o17*o3 + & - 1.78666666666667d3*o11*o15*o16*o17*o3 - 5.6888888888889d1*o14 + & *o17*o18*o3 + 2.97777777777778d2*o15*o17*o18*o3 - 2.133333333333 + & 33d3*o11*o14*o19*o3 + 1.11666666666667d4*o11*o15*o19*o3 - 3.2000 + & 000000000d4*o11*o14*o16*o19*o3 + 1.67500000000000d5*o11*o15*o16* + & o19*o3 + 5.3333333333333d3*o14*o18*o19*o3 - 2.79166666666667d4*o + & 15*o18*o19*o3 - 5.6711111111111d2*o11*o16*o25*o3 - 5.30277777777 + & 78d3*o11*o16*o26*o3 - 3.12500000000000d3*o11*o16*o27*o3 + 2.6666 + & 6666666667d2*o17*o19*o31 + 8.0000000000000d3*o16*o17*o19*o31 - 4 + & .2666666666667d1*o16*o29*o31 + 3.5555555555556d0*o16*o30*o31 - 6 + & .6666666666667d2*o17*o19*o32 + 3.5555555555556d0*o29*o32 - 4.000 + & 0000000000d4*o11*o16*o33 + t3 = 5.0000000000000d3*o18*o33 + 2.50000000000000d5*o16*o18*o33 + & - 6.0000000000000d5*o11*o28*o33 + 1.12500000000000d6*o18*o28*o33 + & - 2.50000000000000d4*o31*o33 - 3.7500000000000d5*o16*o31*o33 + + & 3.12500000000000d4*o32*o33 - 3.9111111111111d1*o10*o16*o18*o4 + + & 1.03111111111111d2*o10*o16*o18*o3*o4 - 4.2000000000000d1*byJ*o14*o + & 6 - 4.0600000000000d2*byJ*o13*o14*o6 + 2.77500000000000d2*byJ*o15*o6 + & + 1.45250000000000d3*byJ*o13*o15*o6 - 3.2375000000000d2*o14*o15*o + & 6 - 4.8241666666667d3*o13*o14*o15*o6 - 2.40000000000000d3*byJ*o16* + & o19*o6 + 2.80000000000000d3*o14*o16*o19*o6 + 2.70666666666667d4* + & o13*o14*o16*o19*o6 - 1.85000000000000d4*o15*o16*o19*o6 - 9.68333 + & 33333333d4*o13*o15*o16*o19*o6 + 1.80000000000000d1*o20*o6 + 2.45 + & 000000000000d1*o21*o6 + 4.7366666666667d2*o13*o21*o6 - 2.3822222 + & 2222222d3*o14*o15*o22*o6 + 2.27555555555556d2*o21*o22*o6 - 1.638 + & 09722222222d4*o14*o15*o23*o6 + 2.28938888888889d3*o21*o23*o6 + t4 = 1.06953125000000d3*o24*o6 + 1.11963541666667d4*o13*o24*o6 + + & 6.2347222222222d3*o22*o24*o6 + 2.93021701388889d4*o23*o24*o6 - + & 1.28000000000000d2*byJ*o14*o3*o6 + 6.7000000000000d2*byJ*o15*o3*o6 - + & 1.76833333333333d3*o14*o15*o3*o6 - 1.27205555555556d4*o13*o14*o + & 15*o3*o6 + 8.5333333333333d3*o14*o16*o19*o3*o6 - 4.4666666666667 + & d4*o15*o16*o19*o3*o6 + 1.49333333333333d2*o21*o3*o6 + 1.44355555 + & 555556d3*o13*o21*o3*o6 + 5.1645833333333d3*o24*o3*o6 + 2.7032638 + & 8888889d4*o13*o24*o3*o6 + 8.0000000000000d4*o28*o33*o6 + 1.02222 + & 222222222d2*o10*o16*o18*o8 - 3.6888888888889d2*o10*o16*o18*o3*o8 + & - 5.6222222222222d2*o11*o16*o4*o8 - 5.3488888888889d3*o11*o16*o + & 22*o4*o8 + 3.5111111111111d3*o11*o16*o3*o4*o8 - 6.6666666666667d + & 1*o10*o16*o18*o9 + 3.3333333333333d2*o10*o16*o18*o3*o9 + 3.66666 + & 66666667d2*o11*o16*o4*o9 + 4.8333333333333d3*o11*o16*o22*o4*o9 - + & 2.80000000000000d3*o11*o16*o3*o4*o9 - 9.5833333333333d2*o11*o16 + & *o8*o9 - 1.72916666666667d4*o11*o16*o22*o8*o9 + 8.2500000000000d + & 3*o11*o16*o3*o8*o9 + ksq(i,j,k) = t1 + t2 + t3 + t4 diff --git a/src/ksq_odd.x b/src/ksq_odd.x new file mode 100644 index 0000000..3f116b9 --- /dev/null +++ b/src/ksq_odd.x @@ -0,0 +1,403 @@ + o1 = mm**2 + o2 = 5.0000000000000d1*o1 + o3 = -8.0000000000000d1 + o2 + o4 = 1/o3 + o5 = 2.00000000000000d0*qgrd(j) + o6 = cos(o5) + o7 = 5.0000000000000d0*o6 + o8 = 3.00000000000000d0 + o7 + o9 = dnnnnngtil + o10 = mm*phigrd(k) + o11 = sin(o10) + o12 = sin(qgrd(j)) + o13 = o12**2 + o14 = o13**2 + o15 = 2.50000000000000d1*o1 + o16 = -4.0000000000000d1 + o15 + o17 = 1/o16 + o18 = dnnngtil + o19 = 2.00000000000000d1*o1 + o20 = -3.2000000000000d1 + o19 + o21 = 1/o20 + o22 = dngtil + o23 = 5.0000000000000d0*o1 + o24 = -8.0000000000000d0 + o23 + o25 = 1/o24 + o26 = cos(o10) + o27 = cos(qgrd(j)) + o28 = 1.03000000000000d2*o1 + o29 = -9.2800000000000d2 + o28 + o30 = o12*o13 + o31 = 3.5000000000000d1*o1 + o32 = -8.3000000000000d1 + o31 + o33 = -1.60000000000000d1 + o1 + o34 = dnnnnnnngtil + o35 = 4.0000000000000d0*qgrd(j) + o36 = cos(o35) + o37 = dnnnngtil + o38 = 1.00000000000000d1*o1 + o39 = -1.60000000000000d1 + o38 + o40 = 1/o39 + o41 = dnngtil + o42 = dnnnnnngtil + o43 = o27**2 + o44 = gtil + o45 = o13*o14 + o46 = 1.07000000000000d2*o1 + o47 = 5.3000000000000d1*o1 + o48 = sin(o5) + o49 = byJ**2 + o50 = o24**2 + o51 = 1/o50 + o52 = o26**2 + o53 = o37**2 + o54 = o1**2 + o55 = o6**2 + o56 = o36**2 + o57 = o39**2 + o58 = 1/o57 + o59 = o41**2 + o60 = o9**2 + o61 = o16**2 + o62 = 1/o61 + o63 = o18**2 + o64 = o20**2 + o65 = 1/o64 + o66 = o22**2 + o67 = o11**2 + o68 = o1*o54 + o69 = o3**2 + o70 = 1/o69 + o71 = o42**2 + o72 = o34**2 + o73 = o14**2 + o74 = o48**2 + o75 = o43**2 + o76 = o44**2 + t1 = -1.00800000000000d2*byJ*o13*o25*o26*o37 - 1.39500000000000d1* + & byJ*o1*o13*o25*o26*o37 - 9.7440000000000d2*byJ*o13*o25*o26*o36*o37 + + & 1.08150000000000d2*byJ*o1*o13*o25*o26*o36*o37 + 1.33200000000000d + & 3*byJ*o13*o26*o40*o41 - 4.6800000000000d2*byJ*o1*o13*o26*o40*o41 + 6 + & .9720000000000d3*byJ*o13*o26*o36*o40*o41 - 2.94000000000000d3*byJ*o1 + & *o13*o26*o36*o40*o41 - 9.6000000000000d2*byJ*o14*o26*o4*o42 + 6.00 + & 00000000000d1*byJ*o1*o14*o26*o4*o42 - 1.80000000000000d2*byJ*o14*o26 + & *o44 + 7.2000000000000d2*byJ*o13*o26*o43*o44 - 2.70000000000000d3* + & byJ*o14*o26*o43*o44 + 4.5000000000000d2*byJ*o26*o44*o45 + 1.80000000 + & 000000d1*o13*o49 - 3.7296000000000d3*o13*o25*o37*o40*o41*o52 + 7 + & .9425000000000d2*o1*o13*o25*o37*o40*o41*o52 - 5.5574400000000d4* + & o13*o25*o36*o37*o40*o41*o52 + 2.21991000000000d4*o1*o13*o25*o36* + & o37*o40*o41*o52 + 2.68800000000000d3*o14*o25*o37*o4*o42*o52 + 2. + & 04000000000000d2*o1*o14*o25*o37*o4*o42*o52 + 2.59840000000000d4* + & o14*o25*o36*o37*o4*o42*o52 - 4.5080000000000d3*o1*o14*o25*o36*o3 + & 7*o4*o42*o52 - 3.5520000000000d4*o14*o4*o40*o41*o42*o52 + 1.4700 + & 0000000000d4*o1*o14*o4*o40*o41*o42*o52 - 1.85920000000000d5*o14* + & o36*o4*o40*o41*o42*o52 + 9.0020000000000d4*o1*o14*o36*o4*o40*o41 + & *o42*o52 + t2 = -1.10400000000000d5*o14*o17*o18*o21*o22*o43*o52 + 2.9160000 + & 0000000d4*o1*o14*o17*o18*o21*o22*o43*o52 + 5.0400000000000d2*o14 + & *o25*o37*o44*o52 + 6.9750000000000d1*o1*o14*o25*o37*o44*o52 + 4. + & 8720000000000d3*o14*o25*o36*o37*o44*o52 - 5.4075000000000d2*o1*o + & 14*o25*o36*o37*o44*o52 - 6.6600000000000d3*o14*o40*o41*o44*o52 + + & 2.34000000000000d3*o1*o14*o40*o41*o44*o52 - 3.4860000000000d4*o + & 14*o36*o40*o41*o44*o52 + 1.47000000000000d4*o1*o14*o36*o40*o41*o + & 44*o52 - 2.01600000000000d3*o13*o25*o37*o43*o44*o52 - 2.79000000 + & 000000d2*o1*o13*o25*o37*o43*o44*o52 + 7.5600000000000d3*o14*o25* + & o37*o43*o44*o52 + 1.04625000000000d3*o1*o14*o25*o37*o43*o44*o52 + & - 1.94880000000000d4*o13*o25*o36*o37*o43*o44*o52 + 2.16300000000 + & 000d3*o1*o13*o25*o36*o37*o43*o44*o52 + 7.3080000000000d4*o14*o25 + & *o36*o37*o43*o44*o52 - 8.1112500000000d3*o1*o14*o25*o36*o37*o43* + & o44*o52 + 2.66400000000000d4*o13*o40*o41*o43*o44*o52 - 9.3600000 + & 000000d3*o1*o13*o40*o41*o43*o44*o52 - 9.9900000000000d4*o14*o40* + & o41*o43*o44*o52 + 3.5100000000000d4*o1*o14*o40*o41*o43*o44*o52 + + & 1.39440000000000d5*o13*o36*o40*o41*o43*o44*o52 - 5.880000000000 + & 0d4*o1*o13*o36*o40*o41*o43*o44*o52 - 5.2290000000000d5*o14*o36*o + & 40*o41*o43*o44*o52 + 2.20500000000000d5*o1*o14*o36*o40*o41*o43*o + & 44*o52 - 1.92000000000000d4*o14*o4*o42*o43*o44*o52 + t3 = 1.20000000000000d3*o1*o14*o4*o42*o43*o44*o52 - 7.6800000000 + & 000d3*o17*o21*o22*o34*o43*o45*o52 + 1.44000000000000d3*o1*o17*o2 + & 1*o22*o34*o43*o45*o52 - 1.26000000000000d3*o25*o37*o44*o45*o52 - + & 1.74375000000000d2*o1*o25*o37*o44*o45*o52 - 1.21800000000000d4* + & o25*o36*o37*o44*o45*o52 + 1.35187500000000d3*o1*o25*o36*o37*o44* + & o45*o52 + 1.66500000000000d4*o40*o41*o44*o45*o52 - 5.85000000000 + & 00d3*o1*o40*o41*o44*o45*o52 + 8.7150000000000d4*o36*o40*o41*o44* + & o45*o52 - 3.6750000000000d4*o1*o36*o40*o41*o44*o45*o52 + 4.80000 + & 00000000d3*o4*o42*o44*o45*o52 - 3.00000000000000d2*o1*o4*o42*o44 + & *o45*o52 + 7.2000000000000d4*o4*o42*o43*o44*o45*o52 - 4.50000000 + & 00000d3*o1*o4*o42*o43*o44*o45*o52 + 1.41120000000000d2*o13*o51*o + & 52*o53 + 3.9060000000000d1*o1*o13*o51*o52*o53 + 2.72832000000000 + & d3*o13*o36*o51*o52*o53 + 7.4760000000000d1*o1*o13*o36*o51*o52*o5 + & 3 + 1.81350000000000d2*o13*o25*o37*o40*o41*o52*o54 - 2.667000000 + & 00000d2*o13*o25*o36*o37*o40*o41*o52*o54 - 2.32500000000000d1*o14 + & *o25*o37*o4*o42*o52*o54 + 1.80250000000000d2*o14*o25*o36*o37*o4* + & o42*o52*o54 - 7.8000000000000d2*o14*o4*o40*o41*o42*o52*o54 - 4.9 + & 000000000000d3*o14*o36*o4*o40*o41*o42*o52*o54 - 1.92000000000000 + & d3*o14*o17*o18*o21*o22*o43*o52*o54 - 6.0000000000000d1*o17*o21*o + & 22*o34*o43*o45*o52*o54 + t4 = 2.70281250000000d0*o13*o51*o52*o53*o54 - 4.1908125000000d1* + & o13*o36*o51*o52*o53*o54 - 2.74432000000000d4*o13*o25*o37*o40*o41 + & *o52*o55 + 9.5432000000000d3*o1*o13*o25*o37*o40*o41*o52*o55 + 1. + & 14688000000000d4*o14*o25*o37*o4*o42*o52*o55 - 4.2560000000000d2* + & o1*o14*o25*o37*o4*o42*o52*o55 - 1.20064000000000d5*o14*o4*o40*o4 + & 1*o42*o52*o55 + 5.2304000000000d4*o1*o14*o4*o40*o41*o42*o52*o55 + & - 1.99200000000000d6*o14*o17*o18*o21*o22*o43*o52*o55 + 2.0850000 + & 0000000d6*o1*o14*o17*o18*o21*o22*o43*o52*o55 + 1.31072000000000d + & 3*o13*o51*o52*o53*o55 + 6.6560000000000d1*o1*o13*o51*o52*o53*o55 + & + 2.60000000000000d2*o13*o25*o37*o40*o41*o52*o54*o55 - 1.820000 + & 00000000d1*o14*o25*o37*o4*o42*o52*o54*o55 - 2.80000000000000d3*o + & 14*o4*o40*o41*o42*o52*o54*o55 - 5.2500000000000d5*o14*o17*o18*o2 + & 1*o22*o43*o52*o54*o55 + 8.4500000000000d-1*o13*o51*o52*o53*o54*o + & 55 - 1.88708800000000d5*o13*o25*o37*o40*o41*o52*o56 + 1.00521050 + & 000000d5*o1*o13*o25*o37*o40*o41*o52*o56 + 1.31868800000000d4*o13 + & *o51*o52*o53*o56 - 2.92726000000000d3*o1*o13*o51*o52*o53*o56 - 8 + & .8322500000000d3*o13*o25*o37*o40*o41*o52*o54*o56 + 1.62450312500 + & 000d2*o13*o51*o52*o53*o54*o56 + 2.46420000000000d4*o13*o52*o58*o + & 59 - 1.73160000000000d4*o1*o13*o52*o58*o59 + 2.57964000000000d5* + & o13*o36*o52*o58*o59 - 1.99416000000000d5*o1*o13*o36*o52*o58*o59 + t5 = 3.04200000000000d3*o13*o52*o54*o58*o59 + 3.8220000000000d4* + & o13*o36*o52*o54*o58*o59 + 1.43648000000000d5*o13*o52*o55*o58*o59 + & - 1.07200000000000d5*o1*o13*o52*o55*o58*o59 + 2.00000000000000d + & 4*o13*o52*o54*o55*o58*o59 + 6.7512200000000d5*o13*o52*o56*o58*o5 + & 9 - 5.6938000000000d5*o1*o13*o52*o56*o58*o59 + 1.20050000000000d + & 5*o13*o52*o54*o56*o58*o59 - 3.07200000000000d2*byJ*o13*o25*o26*o37 + & *o6 - 7.8000000000000d0*byJ*o1*o13*o25*o26*o37*o6 + 3.216000000000 + & 0d3*byJ*o13*o26*o40*o41*o6 - 1.20000000000000d3*byJ*o1*o13*o26*o40*o + & 41*o6 - 1.34400000000000d3*byJ*o14*o26*o4*o42*o6 + 8.4000000000000 + & d1*byJ*o1*o14*o26*o4*o42*o6 - 2.03712000000000d4*o13*o25*o37*o40*o + & 41*o52*o6 + 5.8188000000000d3*o1*o13*o25*o37*o40*o41*o52*o6 - 1. + & 46540800000000d5*o13*o25*o36*o37*o40*o41*o52*o6 + 6.571880000000 + & 0d4*o1*o13*o25*o36*o37*o40*o41*o52*o6 + 1.19552000000000d4*o14*o + & 25*o37*o4*o42*o52*o6 - 1.84000000000000d1*o1*o14*o25*o37*o4*o42* + & o52*o6 + 3.6377600000000d4*o14*o25*o36*o37*o4*o42*o52*o6 - 6.311 + & 2000000000d3*o1*o14*o25*o36*o37*o4*o42*o52*o6 - 1.35488000000000 + & d5*o14*o4*o40*o41*o42*o52*o6 + 5.7940000000000d4*o1*o14*o4*o40*o + & 41*o42*o52*o6 - 2.60288000000000d5*o14*o36*o4*o40*o41*o42*o52*o6 + & + 1.26028000000000d5*o1*o14*o36*o4*o40*o41*o42*o52*o6 + 9.50400 + & 00000000d5*o14*o17*o18*o21*o22*o43*o52*o6 + t6 = -6.3960000000000d5*o1*o14*o17*o18*o21*o22*o43*o52*o6 + 1.53 + & 600000000000d3*o14*o25*o37*o44*o52*o6 + 3.9000000000000d1*o1*o14 + & *o25*o37*o44*o52*o6 - 1.60800000000000d4*o14*o40*o41*o44*o52*o6 + & + 6.0000000000000d3*o1*o14*o40*o41*o44*o52*o6 - 6.1440000000000d + & 3*o13*o25*o37*o43*o44*o52*o6 - 1.56000000000000d2*o1*o13*o25*o37 + & *o43*o44*o52*o6 + 2.30400000000000d4*o14*o25*o37*o43*o44*o52*o6 + & + 5.8500000000000d2*o1*o14*o25*o37*o43*o44*o52*o6 + 6.4320000000 + & 000d4*o13*o40*o41*o43*o44*o52*o6 - 2.40000000000000d4*o1*o13*o40 + & *o41*o43*o44*o52*o6 - 2.41200000000000d5*o14*o40*o41*o43*o44*o52 + & *o6 + 9.0000000000000d4*o1*o14*o40*o41*o43*o44*o52*o6 - 2.688000 + & 00000000d4*o14*o4*o42*o43*o44*o52*o6 + 1.68000000000000d3*o1*o14 + & *o4*o42*o43*o44*o52*o6 + 3.8400000000000d4*o17*o21*o22*o34*o43*o + & 45*o52*o6 - 2.64000000000000d4*o1*o17*o21*o22*o34*o43*o45*o52*o6 + & - 3.8400000000000d3*o25*o37*o44*o45*o52*o6 - 9.7500000000000d1* + & o1*o25*o37*o44*o45*o52*o6 + 4.0200000000000d4*o40*o41*o44*o45*o5 + & 2*o6 - 1.50000000000000d4*o1*o40*o41*o44*o45*o52*o6 + 6.72000000 + & 00000d3*o4*o42*o44*o45*o52*o6 - 4.2000000000000d2*o1*o4*o42*o44* + & o45*o52*o6 + 1.00800000000000d5*o4*o42*o43*o44*o45*o52*o6 - 6.30 + & 00000000000d3*o1*o4*o42*o43*o44*o45*o52*o6 + 8.6016000000000d2*o + & 13*o51*o52*o53*o6 + 1.40880000000000d2*o1*o13*o51*o52*o53*o6 + t7 = 8.3148800000000d3*o13*o36*o51*o52*o53*o6 - 7.1176000000000d + & 2*o1*o13*o36*o51*o52*o53*o6 + 5.6640000000000d2*o13*o25*o37*o40* + & o41*o52*o54*o6 - 2.96800000000000d3*o13*o25*o36*o37*o40*o41*o52* + & o54*o6 - 4.5550000000000d1*o14*o25*o37*o4*o42*o52*o54*o6 + 2.523 + & 50000000000d2*o14*o25*o36*o37*o4*o42*o52*o54*o6 - 3.092000000000 + & 00d3*o14*o4*o40*o41*o42*o52*o54*o6 - 6.8600000000000d3*o14*o36*o + & 4*o40*o41*o42*o52*o54*o6 + 6.9000000000000d4*o14*o17*o18*o21*o22 + & *o43*o52*o54*o6 + 1.50000000000000d3*o17*o21*o22*o34*o43*o45*o52 + & *o54*o6 + 3.02250000000000d0*o13*o51*o52*o53*o54*o6 - 2.34325000 + & 000000d1*o13*o36*o51*o52*o53*o54*o6 + 1.18992000000000d5*o13*o52 + & *o58*o59*o6 - 8.6208000000000d4*o1*o13*o52*o58*o59*o6 + 6.228320 + & 0000000d5*o13*o36*o52*o58*o59*o6 - 4.9504000000000d5*o1*o13*o36* + & o52*o58*o59*o6 + 1.56000000000000d4*o13*o52*o54*o58*o59*o6 + 9.8 + & 000000000000d4*o13*o36*o52*o54*o58*o59*o6 + 6.1952000000000d2*o1 + & 4*o43*o51*o52*o60 - 1.93600000000000d2*o1*o14*o43*o51*o52*o60 + + & 1.51250000000000d1*o14*o43*o51*o52*o54*o60 + 4.3059200000000d3*o + & 14*o43*o51*o52*o55*o60 - 9.5584000000000d2*o1*o14*o43*o51*o52*o5 + & 5*o60 + 5.3045000000000d1*o14*o43*o51*o52*o54*o55*o60 - 3.266560 + & 0000000d3*o14*o43*o51*o52*o6*o60 + 8.7296000000000d2*o1*o14*o43* + & o51*o52*o6*o60 - 5.6650000000000d1*o14*o43*o51*o52*o54*o6*o60 + t8 = 1.47200000000000d4*o18*o34*o43*o45*o52*o62 - 2.968000000000 + & 00d3*o1*o18*o34*o43*o45*o52*o62 + 1.28000000000000d2*o18*o34*o43 + & *o45*o52*o54*o62 - 5.3120000000000d4*o18*o34*o43*o45*o52*o6*o62 + & + 2.57200000000000d4*o1*o18*o34*o43*o45*o52*o6*o62 - 1.400000000 + & 00000d3*o18*o34*o43*o45*o52*o54*o6*o62 + 1.05800000000000d5*o14* + & o43*o52*o62*o63 - 2.94400000000000d4*o1*o14*o43*o52*o62*o63 + 2. + & 04800000000000d3*o14*o43*o52*o54*o62*o63 + 1.37780000000000d6*o1 + & 4*o43*o52*o55*o62*o63 - 1.16200000000000d6*o1*o14*o43*o52*o55*o6 + & 2*o63 + 2.45000000000000d5*o14*o43*o52*o54*o55*o62*o63 - 7.63600 + & 00000000d5*o14*o43*o52*o6*o62*o63 + 4.2824000000000d5*o1*o14*o43 + & *o52*o6*o62*o63 - 4.4800000000000d4*o14*o43*o52*o54*o6*o62*o63 + + & 2.88000000000000d4*o14*o43*o52*o65*o66 - 7.2000000000000d3*o1*o + & 14*o43*o52*o65*o66 + 4.5000000000000d2*o14*o43*o52*o54*o65*o66 + + & 7.2000000000000d5*o14*o43*o52*o55*o65*o66 - 9.0000000000000d5*o + & 1*o14*o43*o52*o55*o65*o66 + 2.81250000000000d5*o14*o43*o52*o54*o + & 55*o65*o66 - 2.88000000000000d5*o14*o43*o52*o6*o65*o66 + 2.16000 + & 000000000d5*o1*o14*o43*o52*o6*o65*o66 - 2.25000000000000d4*o14*o + & 43*o52*o54*o6*o65*o66 - 2.10600000000000d4*o1*o14*o17*o18*o21*o2 + & 2*o67 - 6.7562000000000d4*o1*o14*o17*o18*o22*o25*o67 - 7.0180000 + & 000000d4*o1*o14*o17*o18*o22*o40*o67 + 1.92000000000000d3*o1*o14* + & o17*o42*o43*o44*o67 + t9 = -3.3600000000000d3*o1*o17*o42*o43*o44*o45*o67 - 5.440000000 + & 0000d2*o1*o17*o25*o27*o30*o37*o42*o48*o67 + 1.98400000000000d3*o + & 1*o17*o27*o30*o40*o41*o42*o48*o67 - 1.02000000000000d3*o1*o12*o2 + & 5*o27*o37*o44*o48*o67 + 1.78500000000000d3*o1*o25*o27*o30*o37*o4 + & 4*o48*o67 + 3.7200000000000d3*o1*o12*o27*o40*o41*o44*o48*o67 - 6 + & .5100000000000d3*o1*o27*o30*o40*o41*o44*o48*o67 + 4.625200000000 + & 0d4*o14*o17*o18*o22*o25*o54*o67 + 6.6494000000000d4*o14*o17*o18* + & o22*o40*o54*o67 - 1.20000000000000d2*o14*o17*o42*o43*o44*o54*o67 + & + 2.10000000000000d2*o17*o42*o43*o44*o45*o54*o67 + 1.2200000000 + & 0000d2*o17*o25*o27*o30*o37*o42*o48*o54*o67 - 5.7200000000000d2*o + & 17*o27*o30*o40*o41*o42*o48*o54*o67 + 1.65000000000000d2*o12*o25* + & o27*o37*o44*o48*o54*o67 - 2.88750000000000d2*o25*o27*o30*o37*o44 + & *o48*o54*o67 - 8.4000000000000d2*o12*o27*o40*o41*o44*o48*o54*o67 + & + 1.47000000000000d3*o27*o30*o40*o41*o44*o48*o54*o67 - 5.850000 + & 0000000d4*o1*o14*o17*o18*o21*o22*o55*o67 - 1.42506000000000d5*o1 + & *o14*o17*o18*o22*o25*o55*o67 - 1.31652000000000d5*o1*o14*o17*o18 + & *o22*o40*o55*o67 + 9.2016000000000d4*o14*o17*o18*o22*o25*o54*o55 + & *o67 + 1.25082000000000d5*o14*o17*o18*o22*o40*o54*o55*o67 - 7.02 + & 00000000000d4*o1*o14*o17*o18*o21*o22*o6*o67 - 1.96332000000000d5 + & *o1*o14*o17*o18*o22*o25*o6*o67 - 1.92408000000000d5*o1*o14*o17*o + & 18*o22*o40*o6*o67 - 1.29920000000000d3*o1*o17*o25*o27*o30*o37*o4 + & 2*o48*o6*o67 + 1.33760000000000d4*o1*o17*o27*o30*o40*o41*o42*o48 + & *o6*o67 + t10 = -2.43600000000000d3*o1*o12*o25*o27*o37*o44*o48*o6*o67 + 4. + & 2630000000000d3*o1*o25*o27*o30*o37*o44*o48*o6*o67 + 2.5080000000 + & 0000d4*o1*o12*o27*o40*o41*o44*o48*o6*o67 - 4.3890000000000d4*o1* + & o27*o30*o40*o41*o44*o48*o6*o67 + 1.30692000000000d5*o14*o17*o18* + & o22*o25*o54*o6*o67 + 1.82664000000000d5*o14*o17*o18*o22*o40*o54* + & o6*o67 + 1.08400000000000d2*o17*o25*o27*o30*o37*o42*o48*o54*o6*o + & 67 - 6.5320000000000d3*o17*o27*o30*o40*o41*o42*o48*o54*o6*o67 + + & 5.1000000000000d1*o12*o25*o27*o37*o44*o48*o54*o6*o67 - 8.9250000 + & 000000d1*o25*o27*o30*o37*o44*o48*o54*o6*o67 - 1.06800000000000d4 + & *o12*o27*o40*o41*o44*o48*o54*o6*o67 + 1.86900000000000d4*o27*o30 + & *o40*o41*o44*o48*o54*o6*o67 + 4.0960000000000d3*o1*o14*o60*o62*o + & 67 - 5.1200000000000d2*o14*o54*o60*o62*o67 + 9.2160000000000d3*o + & 1*o14*o55*o60*o62*o67 - 1.15200000000000d3*o14*o54*o55*o60*o62*o + & 67 + 1.22880000000000d4*o1*o14*o6*o60*o62*o67 - 1.53600000000000 + & d3*o14*o54*o6*o60*o62*o67 + 2.63438000000000d5*o1*o14*o62*o63*o6 + & 7 - 1.49158000000000d5*o14*o54*o62*o63*o67 + 5.8030200000000d5*o + & 1*o14*o55*o62*o63*o67 - 3.2531400000000d5*o14*o54*o55*o62*o63*o6 + & 7 + 7.8147600000000d5*o1*o14*o6*o62*o63*o67 - 4.4056800000000d5* + & o14*o54*o6*o62*o63*o67 + 6.8890000000000d3*o1*o14*o51*o66*o67 - + & 5.8100000000000d3*o14*o51*o54*o66*o67 + 1.36890000000000d4*o1*o1 + & 4*o51*o55*o66*o67 + t11 = -1.05300000000000d4*o14*o51*o54*o55*o66*o67 + 1.4641000000 + & 0000d4*o1*o14*o58*o66*o67 - 1.69400000000000d4*o14*o54*o58*o66*o + & 67 + 2.52810000000000d4*o1*o14*o55*o58*o66*o67 - 2.8620000000000 + & 0d4*o14*o54*o55*o58*o66*o67 + 1.94220000000000d4*o1*o14*o51*o6*o + & 66*o67 - 1.56600000000000d4*o14*o51*o54*o6*o66*o67 + 3.847800000 + & 0000d4*o1*o14*o58*o6*o66*o67 - 4.4040000000000d4*o14*o54*o58*o6* + & o66*o67 + 8.1000000000000d3*o1*o14*o65*o66*o67 + 2.2500000000000 + & 0d4*o1*o14*o55*o65*o66*o67 + 2.70000000000000d4*o1*o14*o6*o65*o6 + & 6*o67 - 7.4900000000000d3*o14*o17*o18*o22*o25*o67*o68 - 1.498000 + & 00000000d4*o14*o17*o18*o22*o40*o67*o68 - 5.5000000000000d0*o17*o + & 25*o27*o30*o37*o42*o48*o67*o68 + 2.80000000000000d1*o17*o27*o30* + & o40*o41*o42*o48*o67*o68 - 1.43100000000000d4*o14*o17*o18*o22*o25 + & *o55*o67*o68 - 2.86200000000000d4*o14*o17*o18*o22*o40*o55*o67*o6 + & 8 - 2.07600000000000d4*o14*o17*o18*o22*o25*o6*o67*o68 - 4.152000 + & 0000000d4*o14*o17*o18*o22*o40*o6*o67*o68 - 1.70000000000000d0*o1 + & 7*o25*o27*o30*o37*o42*o48*o6*o67*o68 + 3.5600000000000d2*o17*o27 + & *o30*o40*o41*o42*o48*o6*o67*o68 + 1.60000000000000d1*o14*o60*o62 + & *o67*o68 + 3.6000000000000d1*o14*o55*o60*o62*o67*o68 + 4.8000000 + & 000000d1*o14*o6*o60*o62*o67*o68 + 2.28980000000000d4*o14*o62*o63 + & *o67*o68 + 5.0562000000000d4*o14*o55*o62*o63*o67*o68 + t12 = 6.8052000000000d4*o14*o6*o62*o63*o67*o68 + 1.2250000000000 + & 0d3*o14*o51*o66*o67*o68 + 2.02500000000000d3*o14*o51*o55*o66*o67 + & *o68 + 4.9000000000000d3*o14*o58*o66*o67*o68 + 8.1000000000000d3 + & *o14*o55*o58*o66*o67*o68 + 3.15000000000000d3*o14*o51*o6*o66*o67 + & *o68 + 1.26000000000000d4*o14*o58*o6*o66*o67*o68 + 1.42420000000 + & 000d4*o1*o14*o60*o67*o70 - 1.90400000000000d3*o14*o54*o60*o67*o7 + & 0 + 3.15540000000000d4*o1*o14*o55*o60*o67*o70 - 4.2480000000000d + & 3*o14*o54*o55*o60*o67*o70 + 4.2396000000000d4*o1*o14*o6*o60*o67* + & o70 - 5.6880000000000d3*o14*o54*o6*o60*o67*o70 + 6.4000000000000 + & d1*o14*o60*o67*o68*o70 + 1.44000000000000d2*o14*o55*o60*o67*o68* + & o70 + 1.92000000000000d2*o14*o6*o60*o67*o68*o70 + 5.120000000000 + & 0d2*o1*o43*o45*o62*o67*o71 - 6.4000000000000d1*o43*o45*o54*o62*o + & 67*o71 + 2.00000000000000d0*o43*o45*o62*o67*o68*o71 + 1.28000000 + & 000000d4*o45*o52*o70*o71 - 1.60000000000000d3*o1*o45*o52*o70*o71 + & + 5.0000000000000d1*o45*o52*o54*o70*o71 + 2.50880000000000d4*o4 + & 5*o52*o55*o70*o71 - 3.13600000000000d3*o1*o45*o52*o55*o70*o71 + + & 9.8000000000000d1*o45*o52*o54*o55*o70*o71 + 3.5840000000000d4*o4 + & 5*o52*o6*o70*o71 - 4.4800000000000d3*o1*o45*o52*o6*o70*o71 + 1.4 + & 0000000000000d2*o45*o52*o54*o6*o70*o71 + t13 = -1.20000000000000d4*o4*o42*o44*o52*o73 + 7.5000000000000d2 + & *o1*o4*o42*o44*o52*o73 - 1.68000000000000d4*o4*o42*o44*o52*o6*o7 + & 3 + 1.05000000000000d3*o1*o4*o42*o44*o52*o6*o73 + 5.120000000000 + & 0d2*o43*o52*o62*o72*o73 - 6.4000000000000d1*o1*o43*o52*o62*o72*o + & 73 + 2.00000000000000d0*o43*o52*o54*o62*o72*o73 - 1.054000000000 + & 00d3*o1*o25*o37*o40*o41*o67*o74 + 1.44500000000000d2*o1*o51*o53* + & o67*o74 + 4.0850000000000d2*o25*o37*o40*o41*o54*o67*o74 - 4.6750 + & 000000000d1*o51*o53*o54*o67*o74 - 1.69708000000000d4*o1*o25*o37* + & o40*o41*o55*o67*o74 + 8.2418000000000d2*o1*o51*o53*o55*o67*o74 + + & 7.5821000000000d3*o25*o37*o40*o41*o54*o55*o67*o74 - 3.451000000 + & 0000d1*o51*o53*o54*o55*o67*o74 + 1.92200000000000d3*o1*o58*o59*o + & 67*o74 - 8.6800000000000d2*o54*o58*o59*o67*o74 + 8.7362000000000 + & d4*o1*o55*o58*o59*o67*o74 - 7.4404000000000d4*o54*o55*o58*o59*o6 + & 7*o74 - 9.6232000000000d3*o1*o25*o37*o40*o41*o6*o67*o74 + 6.9020 + & 000000000d2*o1*o51*o53*o6*o67*o74 + 4.7966000000000d3*o25*o37*o4 + & 0*o41*o54*o6*o67*o74 - 1.26100000000000d2*o51*o53*o54*o6*o67*o74 + & + 2.59160000000000d4*o1*o58*o59*o6*o67*o74 - 1.68880000000000d4 + & *o54*o58*o59*o6*o67*o74 - 3.8500000000000d1*o25*o37*o40*o41*o67* + & o68*o74 + 3.7812500000000d0*o51*o53*o67*o68*o74 + t14 = -1.51300000000000d2*o25*o37*o40*o41*o55*o67*o68*o74 + 3.61 + & 25000000000d-1*o51*o53*o55*o67*o68*o74 + 9.8000000000000d1*o58*o + & 59*o67*o68*o74 + 1.58420000000000d4*o55*o58*o59*o67*o68*o74 - 5. + & 0140000000000d2*o25*o37*o40*o41*o6*o67*o68*o74 + 2.3375000000000 + & 0d0*o51*o53*o6*o67*o68*o74 + 2.49200000000000d3*o58*o59*o6*o67*o + & 68*o74 - 3.6000000000000d3*o14*o43*o52*o76 + 4.5000000000000d2*o + & 45*o52*o76 + 2.25000000000000d4*o43*o45*o52*o76 + 1.800000000000 + & 00d3*o1*o13*o43*o67*o76 - 6.3000000000000d3*o1*o14*o43*o67*o76 + + & 5.5125000000000d3*o1*o43*o45*o67*o76 - 2.25000000000000d3*o52*o + & 73*o76 + 2.81250000000000d3*o13*o52*o73*o76 - 3.3750000000000d4* + & o43*o52*o73*o76 + 7.2000000000000d3*o13*o52*o75*o76 - 5.40000000 + & 00000d4*o14*o52*o75*o76 + 1.01250000000000d5*o45*o52*o75*o76 - 1 + & .61920000000000d4*o14*o17*o18*o25*o43*o52*o9 + 4.7828000000000d3 + & *o1*o14*o17*o18*o25*o43*o52*o9 + 8.4480000000000d3*o14*o21*o22*o + & 25*o43*o52*o9 - 2.37600000000000d3*o1*o14*o21*o22*o25*o43*o52*o9 + & - 1.12640000000000d3*o17*o25*o34*o43*o45*o52*o9 + 2.46400000000 + & 000d2*o1*o17*o25*o34*o43*o45*o52*o9 - 3.5200000000000d2*o14*o17* + & o18*o25*o43*o52*o54*o9 + 1.65000000000000d2*o14*o21*o22*o25*o43* + & o52*o54*o9 + t15 = -1.10000000000000d1*o17*o25*o34*o43*o45*o52*o54*o9 - 1.540 + & 48000000000d5*o14*o17*o18*o25*o43*o52*o55*o9 + 8.2058000000000d4 + & *o1*o14*o17*o18*o25*o43*o52*o55*o9 + 1.11360000000000d5*o14*o21* + & o22*o25*o43*o52*o55*o9 - 8.1960000000000d4*o1*o14*o21*o22*o25*o4 + & 3*o52*o55*o9 - 7.2100000000000d3*o14*o17*o18*o25*o43*o52*o54*o55 + & *o9 + 7.7250000000000d3*o14*o21*o22*o25*o43*o52*o54*o55*o9 + 1.0 + & 1120000000000d5*o14*o17*o18*o25*o43*o52*o6*o9 - 4.4447200000000d + & 4*o1*o14*o17*o18*o25*o43*o52*o6*o9 - 6.4512000000000d4*o14*o21*o + & 22*o25*o43*o52*o6*o9 + 3.8256000000000d4*o1*o14*o21*o22*o25*o43* + & o52*o6*o9 + 2.96960000000000d3*o17*o25*o34*o43*o45*o52*o6*o9 - 5 + & .1520000000000d2*o1*o17*o25*o34*o43*o45*o52*o6*o9 + 4.5092000000 + & 000d3*o14*o17*o18*o25*o43*o52*o54*o6*o9 - 4.4340000000000d3*o14* + & o21*o22*o25*o43*o52*o54*o6*o9 + 2.06000000000000d1*o17*o25*o34*o + & 43*o45*o52*o54*o6*o9 + 1.06240000000000d4*o1*o14*o17*o22*o25*o67 + & *o9 - 7.1126000000000d4*o1*o14*o17*o18*o4*o67*o9 + 1.62000000000 + & 000d3*o1*o14*o21*o22*o4*o67*o9 + 2.87980000000000d4*o1*o14*o22*o + & 4*o40*o67*o9 - 5.1440000000000d3*o14*o17*o22*o25*o54*o67*o9 + 3. + & 01060000000000d4*o14*o17*o18*o4*o54*o67*o9 - 1.85960000000000d4* + & o14*o22*o4*o40*o54*o67*o9 + 2.24640000000000d4*o1*o14*o17*o22*o2 + & 5*o55*o67*o9 - 1.52406000000000d5*o1*o14*o17*o18*o4*o55*o67*o9 + + & 4.5000000000000d3*o1*o14*o21*o22*o4*o55*o67*o9 + 5.628600000000 + & 0d4*o1*o14*o22*o4*o40*o55*o67*o9 + t16 = -1.00440000000000d4*o14*o17*o22*o25*o54*o55*o67*o9 + 6.622 + & 2000000000d4*o14*o17*o18*o4*o54*o55*o67*o9 - 3.5676000000000d4*o + & 14*o22*o4*o40*o54*o55*o67*o9 + 3.09120000000000d4*o1*o14*o17*o22 + & *o25*o6*o67*o9 - 2.08212000000000d5*o1*o14*o17*o18*o4*o6*o67*o9 + & + 5.4000000000000d3*o1*o14*o21*o22*o4*o6*o67*o9 + 8.067600000000 + & 0d4*o1*o14*o22*o4*o40*o6*o67*o9 - 1.44120000000000d4*o14*o17*o22 + & *o25*o54*o6*o67*o9 + 8.9304000000000d4*o14*o17*o18*o4*o54*o6*o67 + & *o9 - 5.1648000000000d4*o14*o22*o4*o40*o54*o6*o67*o9 - 5.2096000 + & 000000d4*o1*o14*o18*o62*o67*o9 + 1.69520000000000d4*o14*o18*o54* + & o62*o67*o9 - 1.16928000000000d5*o1*o14*o18*o55*o62*o67*o9 + 3.78 + & 36000000000d4*o14*o18*o54*o55*o62*o67*o9 - 1.56096000000000d5*o1 + & *o14*o18*o6*o62*o67*o9 + 5.0652000000000d4*o14*o18*o54*o6*o62*o6 + & 7*o9 + 2.80000000000000d2*o14*o17*o22*o25*o67*o68*o9 - 1.7120000 + & 0000000d3*o14*o17*o18*o4*o67*o68*o9 + 1.12000000000000d3*o14*o22 + & *o4*o40*o67*o68*o9 + 5.4000000000000d2*o14*o17*o22*o25*o55*o67*o + & 68*o9 - 3.8160000000000d3*o14*o17*o18*o4*o55*o67*o68*o9 + 2.1600 + & 0000000000d3*o14*o22*o4*o40*o55*o67*o68*o9 + 7.8000000000000d2*o + & 14*o17*o22*o25*o6*o67*o68*o9 - 5.1120000000000d3*o14*o17*o18*o4* + & o6*o67*o68*o9 + 3.12000000000000d3*o14*o22*o4*o40*o6*o67*o68*o9 + & - 8.5600000000000d2*o14*o18*o62*o67*o68*o9 - 1.90800000000000d3* + & o14*o18*o55*o62*o67*o68*o9 - 2.55600000000000d3*o14*o18*o6*o62*o + & 67*o68*o9 + ksq(i,j,k) = t1 + t2 + t3 + t4 + t5 + t6 + t7 + t8 + t9 + t10 + + & t11 + t12 + t13 + t14 + t15 + t16 diff --git a/src/make.code.defn b/src/make.code.defn new file mode 100644 index 0000000..da14e77 --- /dev/null +++ b/src/make.code.defn @@ -0,0 +1,9 @@ +# Main make.code.defn file for thorn RotatingDBHIVP +# $Header$ + +# Source files in this directory +SRCS = RotatingDBHIVP.F Stab3d.F + +# Subdirectories containing source files +SUBDIRS = + diff --git a/src/psi_1st_deriv.x b/src/psi_1st_deriv.x new file mode 100644 index 0000000..c7725c5 --- /dev/null +++ b/src/psi_1st_deriv.x @@ -0,0 +1,20 @@ + o1 = 5.0000000000000d-1*eta(i,j,k) + o2 = exp(o1) + o3 = psi3d(i,j,k) + o4 = 1/o3 + o5 = cos(phi(i,j,k)) + o6 = cos(q(i,j,k)) + o7 = dqpsi3d(i,j,k) + o8 = -1.50000000000000d0*eta(i,j,k) + o9 = exp(o8) + o10 = dphipsi3d(i,j,k) + o11 = sin(phi(i,j,k)) + o12 = sin(q(i,j,K)) + o13 = 1/o12 + o14 = detapsi3d(i,j,k) + psix(i,j,k) = o2*o4*(-(o10*o11*o13*o9) + o12*o14*o5*o9 - 5.00000 + & 00000000d-1*o12*o3*o5*o9 + o5*o6*o7*o9) + psiy(i,j,k) = o2*o4*(o11*o12*o14*o9 - 5.0000000000000d-1*o11*o12 + & *o3*o9 + 1.00000000000000d0*o10*o13*o5*o9 + o11*o6*o7*o9) + psiz(i,j,k) = o2*o4*(o14*o6*o9 - 5.0000000000000d-1*o3*o6*o9 - o + & 12*o7*o9) diff --git a/src/psi_2nd_deriv.x b/src/psi_2nd_deriv.x new file mode 100644 index 0000000..8cc3e60 --- /dev/null +++ b/src/psi_2nd_deriv.x @@ -0,0 +1,71 @@ + o1 = 5.0000000000000d-1*eta(i,j,k) + o2 = exp(o1) + o3 = psi3d(i,j,k) + o4 = 1/o3 + o5 = cos(phi(i,j,k)) + o6 = o5**2 + o7 = cos(q(i,j,k)) + o8 = o7**2 + o9 = detapsi3d(i,j,k) + o10 = -2.50000000000000d0*eta(i,j,k) + o11 = exp(o10) + o12 = dqqpsi3d(i,j,k) + o13 = detaphipsi3d(i,j,k) + o14 = sin(phi(i,j,k)) + o15 = dphipsi3d(i,j,k) + o16 = o14**2 + o17 = sin(q(i,j,k)) + o18 = o17**2 + o19 = 1/o18 + o20 = dphiphipsi3d(i,j,k) + o21 = detaqpsi3d(i,j,k) + o22 = dqpsi3d(i,j,k) + o23 = detaetapsi3d(i,j,k) + o24 = tan(q(i,j,k)) + o25 = o24**2 + o26 = 1/o25 + o27 = dqphipsi3d(i,j,k) + o28 = 1/o24 + psixx(i,j,k) = o2*o4*(1.00000000000000d0*o11*o16*o19*o20 + 1.000 + & 00000000000d0*o11*o16*o22*o28 - 5.0000000000000d-1*o11*o16*o3 - + & 2.00000000000000d0*o11*o13*o14*o5 + 2.00000000000000d0*o11*o14*o + & 15*o5 + 1.00000000000000d0*o11*o14*o15*o19*o5 + 1.00000000000000 + & d0*o11*o14*o15*o26*o5 - 2.00000000000000d0*o11*o14*o27*o28*o5 + + & o11*o18*o23*o6 + 7.5000000000000d-1*o11*o18*o3*o6 + 2.0000000000 + & 0000d0*o11*o17*o21*o6*o7 - 3.00000000000000d0*o11*o17*o22*o6*o7 + & + o11*o12*o6*o8 - 5.0000000000000d-1*o11*o3*o6*o8 + o11*o16*o9 - + & 2.00000000000000d0*o11*o18*o6*o9 + o11*o6*o8*o9) + psixy(i,j,k) = o2*o4*(-(o11*o13*o16) + 1.50000000000000d0*o11*o1 + & 5*o16 + 1.00000000000000d0*o11*o15*o16*o26 - o11*o16*o27*o28 - o + & 11*o14*o19*o20*o5 + o11*o14*o18*o23*o5 - o11*o14*o22*o28*o5 + 5. + & 0000000000000d-1*o11*o14*o3*o5 + 7.5000000000000d-1*o11*o14*o18* + & o3*o5 + o11*o13*o6 - 5.0000000000000d-1*o11*o15*o6 - o11*o15*o19 + & *o6 + 1.00000000000000d0*o11*o27*o28*o6 + 2.00000000000000d0*o11 + & *o14*o17*o21*o5*o7 - 3.00000000000000d0*o11*o14*o17*o22*o5*o7 + + & o11*o12*o14*o5*o8 - 5.0000000000000d-1*o11*o14*o3*o5*o8 - o11*o1 + & 4*o5*o9 - 2.00000000000000d0*o11*o14*o18*o5*o9 + o11*o14*o5*o8*o + & 9) + psixz(i,j,k) = o2*o4*(o11*o14*o27 - o11*o13*o14*o28 + 5.00000000 + & 00000d-1*o11*o14*o15*o28 - o11*o18*o21*o5 + 1.50000000000000d0*o + & 11*o18*o22*o5 - o11*o12*o17*o5*o7 + o11*o17*o23*o5*o7 + 1.250000 + & 00000000d0*o11*o17*o3*o5*o7 + o11*o21*o5*o8 - 1.50000000000000d0 + & *o11*o22*o5*o8 - 3.00000000000000d0*o11*o17*o5*o7*o9) + psiyy(i,j,k) = o2*o4*(o11*o16*o18*o23 + 7.5000000000000d-1*o11*o + & 16*o18*o3 + 2.00000000000000d0*o11*o13*o14*o5 - 2.00000000000000 + & d0*o11*o14*o15*o5 - o11*o14*o15*o19*o5 - o11*o14*o15*o26*o5 + 2. + & 00000000000000d0*o11*o14*o27*o28*o5 + 1.00000000000000d0*o11*o19 + & *o20*o6 + 1.00000000000000d0*o11*o22*o28*o6 - 5.0000000000000d-1 + & *o11*o3*o6 + 2.00000000000000d0*o11*o16*o17*o21*o7 - 3.000000000 + & 00000d0*o11*o16*o17*o22*o7 + o11*o12*o16*o8 - 5.0000000000000d-1 + & *o11*o16*o3*o8 - 2.00000000000000d0*o11*o16*o18*o9 + o11*o6*o9 + + & o11*o16*o8*o9) + psiyz(i,j,k) = o2*o4*(-(o11*o14*o18*o21) + 1.50000000000000d0*o1 + & 1*o14*o18*o22 - o11*o27*o5 + 1.00000000000000d0*o11*o13*o28*o5 - + & 5.0000000000000d-1*o11*o15*o28*o5 - o11*o12*o14*o17*o7 + o11*o1 + & 4*o17*o23*o7 + 1.25000000000000d0*o11*o14*o17*o3*o7 + o11*o14*o2 + & 1*o8 - 1.50000000000000d0*o11*o14*o22*o8 - 3.00000000000000d0*o1 + & 1*o14*o17*o7*o9) + psizz(i,j,k) = o2*o4*(o11*o12*o18 - 5.0000000000000d-1*o11*o18*o + & 3 - 2.00000000000000d0*o11*o17*o21*o7 + 3.00000000000000d0*o11*o + & 17*o22*o7 + o11*o23*o8 + 7.5000000000000d-1*o11*o3*o8 + o11*o18* + & o9 - 2.00000000000000d0*o11*o8*o9) |