Menu
▾
▴
[Maxima-commits] CVS: maxima/share/contrib/gentran/test README,NONE,1.1 array.mac,NONE,1.1 array.output,NONE,1.1 datatest.mac,NONE,1.1 datatest.out,NONE,1.1 for.mac,NONE,1.1 graeffe.f,NONE,1.1 graeffe.mac,NONE,1.1 graeffe.output,NONE,1.1 ham.in,NONE,1.1 loop.mac,NONE,1.1 matrix.f,NONE,1.1 matrix.mac,NONE,1.1 matrix.output,NONE,1.1 mcond.lisp,NONE,1.1 mdo.lisp,NONE,1.1 rk.in,NONE,1.1 runge.mac,NONE,1.1 runge.output,NONE,1.1 runge.template,NONE,1.1 t.mac,NONE,1.1 type.mac,NONE,1.1 type.output,NONE,1.1
Update of /cvsroot/maxima/maxima/share/contrib/gentran/test
In directory sc8-pr-cvs1:/tmp/cvs-serv31892/test
Added Files:
README array.mac array.output datatest.mac datatest.out
for.mac graeffe.f graeffe.mac graeffe.output ham.in loop.mac
matrix.f matrix.mac matrix.output mcond.lisp mdo.lisp rk.in
runge.mac runge.output runge.template t.mac type.mac
type.output
Log Message:
Initial import of gentran
--- NEW FILE: README ---
Files here can be used to test the gentran package. It include some
example programs described in section 6 of the the gentran manual.
The name.v file can be batched into vaxima and a name.f file will be produced.
This should match exactly the corresponding file name.output.
To perform the test:
invoke vaxima
(C1) loadfile("../gtload.l")$
(C2) batch("filename.v")$
There are some .l files that can be loaded into the lisp level
under Vaxima to test the lisp level calls to gentran.
--- NEW FILE: array.mac ---
mat : genmatrix(mat, 3,3, 1,1)$
mat[1,1] :
18*cos(q3)*cos(q2)*m30*p^2 - 9*sin(q3)^2*p^2*m30
- sin(q3)^2*j30y + sin(q3)^2*j30z + p^2*m10 + 18*p^2*m30
+ j10y + j30y$
mat[2,1] :
mat[1,2] :
9*cos(q3)*cos(q2)*m30*p^2 - sin(q3)^2*j30y + sin(q3)^2*j30z
- 9*sin(q3)^2*m30*p^2 + j30y + 9*m30*p^2$
mat[3,1] :
mat[1,3] :
-9*sin(q3)*sin(q2)*m30*p^2$
mat[2,2] :
-sin(q3)^2*j30y + sin(q3)^2*j30z - 9*sin(q3)^2 *m30*p^2
+ j30y + 9*m30*p^2$
mat[3,2] : mat[2,3] : 0$
mat[3,3] : 9*m30*p^2 + j30x$
gentranout("array.f")$
for i:1 thru 3 do
for j:i thru 3 do
gentran( lrsetq(mat[i,j], mat[i,j]) )$
gentranshut("array.f")$
--- NEW FILE: array.output ---
mat(1,1)=j10y+j30y+m10*p**2+18.0*m30*p**2+18.0*m30*p**2*cos(q2)*
. cos(q3)-j30y*sin(q3)**2+j30z*sin(q3)**2-9.0*m30*p**2*sin(q3)**2
mat(1,2)=j30y+9.0*m30*p**2+9.0*m30*p**2*cos(q2)*cos(q3)-j30y*sin(
. q3)**2+j30z*sin(q3)**2-9.0*m30*p**2*sin(q3)**2
mat(1,3)=-9.0*m30*p**2*sin(q2)*sin(q3)
mat(2,2)=j30y+9.0*m30*p**2-j30y*sin(q3)**2+j30z*sin(q3)**2-9.0*m30
. *p**2*sin(q3)**2
mat(2,3)=0
mat(3,3)=j30x+9.0*m30*p**2
--- NEW FILE: datatest.mac ---
gentran(data("a(10)", 3.1416159, 9.0, 8.1728907, 0.009988, 0.99887766, 8.7906, 6.5432, 1.2345, 0.12345, 98.098, 7654.0987, 11.1111111));
--- NEW FILE: datatest.out ---
(c12) gentran(data("a(10)", 3.1416159, 9.0, 8.1728907, 0.009988, 0.99887766,
8.7906, 6.5432, 1.2345, 0.12345, 98.098, 7654.0987, 11.1111111))
data a(10)/3.1416159,9.0,8.1728907,0.009988,0.99887766,8.7906,
. 6.5432,1.2345,0.12345,98.098,7654.0987,11.1111111/
--- NEW FILE: for.mac ---
gentran(
for i:1 thru 10 do
for j:1 thru 20 do
if i=j
then m[i,j] : 1.0
else m[i,j] : 0.0);
--- NEW FILE: graeffe.f ---
subroutine graeff(a,b)
real*8 a(11),b(11)
c
c -- Graeffe Root-Squaring Method --
c -- to Find Roots of a Polynomial --
c
b(1)=a(11)**2
b(2)=-a(10)**2+2.0*a(9)*a(11)
b(3)=a(9)**2-2.0*a(8)*a(10)+2.0*a(7)*a(11)
b(4)=-a(8)**2+2.0*a(7)*a(9)-2.0*a(6)*a(10)+2.0*a(5)*a(11)
b(5)=a(7)**2-2.0*a(6)*a(8)+2.0*a(5)*a(9)-2.0*a(4)*a(10)+2.0*a(3)*a
. (11)
b(6)=-a(6)**2+2.0*a(5)*a(7)-2.0*a(4)*a(8)+2.0*a(3)*a(9)-2.0*a(2)*a
. (10)+2.0*a(1)*a(11)
b(7)=a(5)**2-2.0*a(4)*a(6)+2.0*a(3)*a(7)-2.0*a(2)*a(8)+2.0*a(1)*a(
. 9)
b(8)=-a(4)**2+2.0*a(3)*a(5)-2.0*a(2)*a(6)+2.0*a(1)*a(7)
b(9)=a(3)**2-2.0*a(2)*a(4)+2.0*a(1)*a(5)
b(10)=-a(2)**2+2.0*a(1)*a(3)
b(11)=a(1)**2
return
end
--- NEW FILE: graeffe.mac ---
n : 10$
q : 0$
for i:0 step 2 thru n do
q : q + a(i+1)*x^(n-i)$
r : 0$
for i:1 step 2 thru n-1 do
r : r + a(i+1)*x^(n-i)$
p : q^2 - r^2$
p : ratsubst(y, x^2, p)$
for i:0 thru n do
b[i] : ratcoeff(p, y, i)$
gentranlang : fortran$
gentranpush("graeffe.f");
gentran( subroutine( graeff(a,b) ),
type("real*8", a(eval(n+1)), b(eval(n+1))),
literal("c",cr),
literal("c",tab,"-- Graeffe Root-Squaring Method --",cr),
literal("c",tab,"-- to Find Roots of a Polynomial --",cr),
literal("c",cr) );
for i:1 thru n+1 do
gentran( lrsetq( b[i], b[i-1] ) );
gentran( return(), end() );
gentranpop(false);
--- NEW FILE: graeffe.output ---
subroutine graeff(a,b)
real*8 a(11),b(11)
c
c -- Graeffe Root-Squaring Method --
c -- to Find Roots of a Polynomial --
c
b(1)=a(11)**2
b(2)=-a(10)**2+2.0*a(9)*a(11)
b(3)=a(9)**2-2.0*a(8)*a(10)+2.0*a(7)*a(11)
b(4)=-a(8)**2+2.0*a(7)*a(9)-2.0*a(6)*a(10)+2.0*a(5)*a(11)
b(5)=a(7)**2-2.0*a(6)*a(8)+2.0*a(5)*a(9)-2.0*a(4)*a(10)+2.0*a(3)*a
. (11)
b(6)=-a(6)**2+2.0*a(5)*a(7)-2.0*a(4)*a(8)+2.0*a(3)*a(9)-2.0*a(2)*a
. (10)+2.0*a(1)*a(11)
b(7)=a(5)**2-2.0*a(4)*a(6)+2.0*a(3)*a(7)-2.0*a(2)*a(8)+2.0*a(1)*a(
. 9)
b(8)=-a(4)**2+2.0*a(3)*a(5)-2.0*a(2)*a(6)+2.0*a(1)*a(7)
b(9)=a(3)**2-2.0*a(2)*a(4)+2.0*a(1)*a(5)
b(10)=-a(2)**2+2.0*a(1)*a(3)
b(11)=a(1)**2
return
end
--- NEW FILE: ham.in ---
/* Hamiltonian Calculation */
difq : diff(h, p)$
difp : -diff(h, q) - ratsubst(p/m, qdot, diff(d, qdot))$
rungekutta(difp, difq, p, q, tt)$
--- NEW FILE: loop.mac ---
gentranlang:fortran;
gentran( type("implicit real*8", "a-h","o-z"),
type("real*8", m(4,4)),
for i:1 thru 4 do
for j:1 thru 4 do
if i=j
then m[i,j] : 1.0
else m[i,j] : 0.0,
type("integer", i,j));
--- NEW FILE: matrix.f ---
c
c --- Calculate Matrix Values ---
c
mat(1,1)=j10y+j30y+m10*p**2+18.0*m30*p**2+18.0*m30*p**2*cos(q2)*
. cos(q3)-j30y*sin(q3)**2+j30z*sin(q3)**2-9.0*m30*p**2*sin(q3)**2
mat(1,2)=j30y+9.0*m30*p**2+9.0*m30*p**2*cos(q2)*cos(q3)-j30y*sin(
. q3)**2+j30z*sin(q3)**2-9.0*m30*p**2*sin(q3)**2
mat(1,3)=-9.0*m30*p**2*sin(q2)*sin(q3)
mat(2,2)=j30y+9.0*m30*p**2-j30y*sin(q3)**2+j30z*sin(q3)**2-9.0*m30
. *p**2*sin(q3)**2
mat(2,3)=0
mat(3,3)=j30x+9.0*m30*p**2
--- NEW FILE: matrix.mac ---
mat : genmatrix(mat, 3,3, 1,1)$
mat[1,1] :
18*cos(q3)*cos(q2)*m30*p^2 - 9*sin(q3)^2*p^2*m30
- sin(q3)^2*j30y + sin(q3)^2*j30z + p^2*m10 + 18*p^2*m30
+ j10y + j30y$
mat[2,1] :
mat[1,2] :
9*cos(q3)*cos(q2)*m30*p^2 - sin(q3)^2*j30y + sin(q3)^2*j30z
- 9*sin(q3)^2*m30*p^2 + j30y + 9*m30*p^2$
mat[3,1] :
mat[1,3] :
-9*sin(q3)*sin(q2)*m30*p^2$
mat[2,2] :
-sin(q3)^2*j30y + sin(q3)^2*j30z - 9*sin(q3)^2 *m30*p^2
+ j30y + 9*m30*p^2$
mat[3,2] : mat[2,3] : 0$
mat[3,3] : 9*m30*p^2 + j30x$
?maxexpprintlen\* : 300$ /* expression segmentation length */
gentranout("matrix.f")$
gentran( literal("c", cr,
"c --- Calculate Matrix Values ---", cr,
"c", cr) )$
for i:1 thru 3 do
for j:i thru 3 do
gentran( lrsetq(mat[i,j], mat[i,j]) )$
--- NEW FILE: matrix.output ---
c
c --- Calculate Matrix Values ---
c
t0=p**2
t1=sin(q3)**2
mat(1,1)=j10y+j30y+m10*t0+18.0*m30*t0+18.0*m30*t0*cos(q2)*cos(q3)-
. j30y*t1+j30z*t1-9.0*m30*t0*t1
t0=p**2
t1=sin(q3)**2
mat(1,2)=j30y+9.0*m30*t0+9.0*m30*t0*cos(q2)*cos(q3)-j30y*t1+j30z*
. t1-9.0*m30*t0*t1
mat(1,3)=-9.0*m30*p**2*sin(q2)*sin(q3)
t0=p**2
t1=sin(q3)**2
mat(2,2)=j30y+9.0*m30*t0-j30y*t1+j30z*t1-9.0*m30*t0*t1
mat(2,3)=0.0
mat(3,3)=j30x+9.0*m30*p**2
c
c --- Calculate Inverse Matrix Values ---
c
t0=-j30x*j30y
t1=-9.0*j30x
t2=-9.0*j30y
t3=p**2
t4=m30**2
t5=p**4
t6=9.0*j30x
t7=9.0*j30y
t8=sin(q3)
t9=t8**2
t10=9.0*j10y+t6
t11=81.0*j30y
t12=m30**3
t13=p**6
t14=-729.0*t12*t13
t15=sin(q2)**2
tt0=j10y*j30x*j30y+(j30x*j30y*m10+(9.0*j10y*j30x+t10*j30y)*m30)*t3
. +((t6+t7)*m10*m30+(81.0*j10y+81.0*j30x+t11)*t4)*t5+(81.0*m10*t4+
. 729.0*t12)*t13+(-81.0*j30x*t4*t5+t14)*cos(q2)**2*cos(q3)**2
tt0=tt0+(-j10y*j30x*j30y+j10y*j30x*j30z+((t0+j30x*j30z)*m10+(-9.0*
. j10y*j30x+(-9.0*j10y+t1)*j30y+t10*j30z)*m30)*t3+((t1+t2+9.0*j30z)
. *m10*m30+(-81.0*j10y-81.0*j30x-81.0*j30y+81.0*j30z)*t4)*t5+(-81.0
. *m10*t4-729.0*t12)*t13+(-81.0*j30y*t4*t5+t14)*t15)*t9
inv(1,1)=-(t0+(t1+t2)*m30*t3-81.0*t4*t5+(j30x*j30y-j30x*j30z+(t6+
. t7-9.0*j30z)*m30*t3+81.0*t4*t5)*t9)/(tt0+((t11-81.0*j30z)*t4*t5+
. 729.0*t12*t13)*t15*t8**4)
t0=-j30x*j30y
t1=-9.0*j30x
t2=-9.0*j30y
t3=p**2
t4=m30**2
t5=p**4
t6=-81.0*t4*t5
t7=cos(q2)
t8=cos(q3)
t9=9.0*j30x
t10=9.0*j30y
t11=sin(q3)
t12=t11**2
t13=9.0*j10y+t9
t14=81.0*j30y
t15=m30**3
t16=p**6
t17=-729.0*t15*t16
t18=sin(q2)**2
tt0=j10y*j30x*j30y+(j30x*j30y*m10+(9.0*j10y*j30x+t13*j30y)*m30)*t3
. +((t9+t10)*m10*m30+(81.0*j10y+81.0*j30x+t14)*t4)*t5+(81.0*m10*t4+
. 729.0*t15)*t16+(-81.0*j30x*t4*t5+t17)*t7**2*t8**2
tt0=tt0+(-j10y*j30x*j30y+j10y*j30x*j30z+((t0+j30x*j30z)*m10+(-9.0*
. j10y*j30x+(-9.0*j10y+t1)*j30y+t13*j30z)*m30)*t3+((t1+t2+9.0*j30z)
. *m10*m30+(-81.0*j10y-81.0*j30x-81.0*j30y+81.0*j30z)*t4)*t5+(-81.0
. *m10*t4-729.0*t15)*t16+(-81.0*j30y*t4*t5+t17)*t18)*t12
inv(1,2)=(t0+(t1+t2)*m30*t3+t6+(-9.0*j30x*m30*t3+t6)*t7*t8+(j30x*
. j30y-j30x*j30z+(t9+t10-9.0*j30z)*m30*t3+81.0*t4*t5)*t12)/(tt0+((
. t14-81.0*j30z)*t4*t5+729.0*t15*t16)*t18*t11**4)
t0=p**2
t1=m30**2
t2=p**4
t3=sin(q2)
t4=sin(q3)
t5=9.0*j30y
t6=9.0*j30x
t7=9.0*j10y+t6
t8=81.0*j30y
t9=m30**3
t10=p**6
t11=-729.0*t9*t10
t12=-9.0*j30x
t13=t3**2
tt0=-j10y*j30x*j30y+j10y*j30x*j30z+((-j30x*j30y+j30x*j30z)*m10+(
. -9.0*j10y*j30x+(-9.0*j10y+t12)*j30y+t7*j30z)*m30)*t0+((t12-9.0*
. j30y+9.0*j30z)*m10*m30+(-81.0*j10y-81.0*j30x-81.0*j30y+81.0*j30z)
. *t1)*t2+(-81.0*m10*t1-729.0*t9)*t10
tt0=j10y*j30x*j30y+(j30x*j30y*m10+(9.0*j10y*j30x+t7*j30y)*m30)*t0+
. ((t6+t5)*m10*m30+(81.0*j10y+81.0*j30x+t8)*t1)*t2+(81.0*m10*t1+
. 729.0*t9)*t10+(-81.0*j30x*t1*t2+t11)*cos(q2)**2*cos(q3)**2+(tt0+(
. -81.0*j30y*t1*t2+t11)*t13)*t4**2
inv(1,3)=-((-9.0*j30y*m30*t0-81.0*t1*t2)*t3*t4+((t5-9.0*j30z)*m30*
. t0+81.0*t1*t2)*t3*t4**3)/(tt0+((t8-81.0*j30z)*t1*t2+729.0*t9*t10)
. *t13*t4**4)
t0=-j30x*j30y
t1=-9.0*j10y
t2=-9.0*j30y
t3=p**2
t4=m30**2
t5=p**4
t6=cos(q2)
t7=cos(q3)
t8=9.0*j30x
t9=9.0*j30y
t10=sin(q2)**2
t11=sin(q3)
t12=t11**2
t13=9.0*j10y+t8
t14=81.0*j30y
t15=m30**3
t16=p**6
t17=-729.0*t15*t16
t18=-9.0*j30x
tt1=-(-j10y*j30x+t0+(-j30x*m10+(t1-18.0*j30x+t2)*m30)*t3+(-9.0*m10
. *m30-162.0*t4)*t5+(-18.0*j30x*m30*t3-162.0*t4*t5)*t6*t7+(j30x*
. j30y-j30x*j30z+(t8+t9-9.0*j30z)*m30*t3+81.0*t4*t5+81.0*t4*t5*t10)
. *t12)
tt0=j10y*j30x*j30y+(j30x*j30y*m10+(9.0*j10y*j30x+t13*j30y)*m30)*t3
. +((t8+t9)*m10*m30+(81.0*j10y+81.0*j30x+t14)*t4)*t5+(81.0*m10*t4+
. 729.0*t15)*t16+(-81.0*j30x*t4*t5+t17)*t6**2*t7**2
tt0=tt0+(-j10y*j30x*j30y+j10y*j30x*j30z+((t0+j30x*j30z)*m10+(-9.0*
. j10y*j30x+(t1+t18)*j30y+t13*j30z)*m30)*t3+((t18+t2+9.0*j30z)*m10*
. m30+(-81.0*j10y-81.0*j30x-81.0*j30y+81.0*j30z)*t4)*t5+(-81.0*m10*
. t4-729.0*t15)*t16+(-81.0*j30y*t4*t5+t17)*t10)*t12
inv(2,2)=tt1/(tt0+((t14-81.0*j30z)*t4*t5+729.0*t15*t16)*t10*t11**4
. )
t0=p**2
t1=m30**2
t2=p**4
t3=sin(q2)
t4=cos(q2)
t5=cos(q3)
t6=sin(q3)
t7=9.0*j30y
t8=9.0*j30x
t9=9.0*j10y+t8
t10=81.0*j30y
t11=m30**3
t12=p**6
t13=-729.0*t11*t12
t14=-9.0*j30x
t15=t3**2
tt0=-j10y*j30x*j30y+j10y*j30x*j30z+((-j30x*j30y+j30x*j30z)*m10+(
. -9.0*j10y*j30x+(-9.0*j10y+t14)*j30y+t9*j30z)*m30)*t0+((t14-9.0*
. j30y+9.0*j30z)*m10*m30+(-81.0*j10y-81.0*j30x-81.0*j30y+81.0*j30z)
. *t1)*t2+(-81.0*m10*t1-729.0*t11)*t12
tt0=j10y*j30x*j30y+(j30x*j30y*m10+(9.0*j10y*j30x+t9*j30y)*m30)*t0+
. ((t8+t7)*m10*m30+(81.0*j10y+81.0*j30x+t10)*t1)*t2+(81.0*m10*t1+
. 729.0*t11)*t12+(-81.0*j30x*t1*t2+t13)*t4**2*t5**2+(tt0+(-81.0*
. j30y*t1*t2+t13)*t15)*t6**2
inv(2,3)=(((-9.0*j30y*m30*t0-81.0*t1*t2)*t3-81.0*t1*t2*t4*t3*t5)*
. t6+((t7-9.0*j30z)*m30*t0+81.0*t1*t2)*t3*t6**3)/(tt0+((t10-81.0*
. j30z)*t1*t2+729.0*t11*t12)*t15*t6**4)
t0=-9.0*j10y
t1=-9.0*j30y
t2=p**2
t3=m30**2
t4=p**4
t5=cos(q2)**2
t6=cos(q3)**2
t7=9.0*j10y
t8=9.0*j30y
t9=sin(q3)
t10=t9**2
t11=9.0*j30x
t12=t7+t11
t13=81.0*j30y
t14=m30**3
t15=p**6
t16=-729.0*t14*t15
t17=-9.0*j30x
t18=sin(q2)**2
tt0=j10y*j30x*j30y+(j30x*j30y*m10+(9.0*j10y*j30x+t12*j30y)*m30)*t2
. +((t11+t8)*m10*m30+(81.0*j10y+81.0*j30x+t13)*t3)*t4+(81.0*m10*t3+
. 729.0*t14)*t15+(-81.0*j30x*t3*t4+t16)*t5*t6
tt0=tt0+(-j10y*j30x*j30y+j10y*j30x*j30z+((-j30x*j30y+j30x*j30z)*
. m10+(-9.0*j10y*j30x+(t0+t17)*j30y+t12*j30z)*m30)*t2+((t17+t1+9.0*
. j30z)*m10*m30+(-81.0*j10y-81.0*j30x-81.0*j30y+81.0*j30z)*t3)*t4+(
. -81.0*m10*t3-729.0*t14)*t15+(-81.0*j30y*t3*t4+t16)*t18)*t10
inv(3,3)=-(-j10y*j30y+(-j30y*m10+(t0+t1)*m30)*t2+(-9.0*m10*m30
. -81.0*t3)*t4+81.0*t3*t4*t5*t6+(j10y*j30y-j10y*j30z+((j30y-j30z)*
. m10+(t7+t8-9.0*j30z)*m30)*t2+(9.0*m10*m30+81.0*t3)*t4)*t10)/(tt0+
. ((t13-81.0*j30z)*t3*t4+729.0*t14*t15)*t18*t9**4)
c
c --- Copy Entries Across Main Diagonals ---
c
do 25001 i=1,3
do 25002 j=i+1.0,3
mat(j,i)=mat(i,j)
inv(j,i)=inv(i,j)
25002 continue
25001 continue
c
c --- Calculate Matrix Values ---
c
t0=p**2
t1=sin(q3)**2
mat(1,1)=j10y+j30y+m10*t0+18.0*m30*t0+18.0*m30*t0*cos(q2)*cos(q3)-
. j30y*t1+j30z*t1-9.0*m30*t0*t1
t0=p**2
t1=sin(q3)**2
mat(1,2)=j30y+9.0*m30*t0+9.0*m30*t0*cos(q2)*cos(q3)-j30y*t1+j30z*
. t1-9.0*m30*t0*t1
mat(1,3)=-9.0*m30*p**2*sin(q2)*sin(q3)
t0=p**2
t1=sin(q3)**2
mat(2,2)=j30y+9.0*m30*t0-j30y*t1+j30z*t1-9.0*m30*t0*t1
mat(2,3)=0.0
mat(3,3)=j30x+9.0*m30*p**2
c
c --- Calculate Inverse Matrix Values ---
c
t0=-j30x*j30y
t1=-9.0*j30x
t2=-9.0*j30y
t3=p**2
t4=m30**2
t5=p**4
t6=9.0*j30x
t7=9.0*j30y
t8=sin(q3)
t9=t8**2
t10=9.0*j10y+t6
t11=81.0*j30y
t12=m30**3
t13=p**6
t14=-729.0*t12*t13
t15=sin(q2)**2
tt0=j10y*j30x*j30y+(j30x*j30y*m10+(9.0*j10y*j30x+t10*j30y)*m30)*t3
. +((t6+t7)*m10*m30+(81.0*j10y+81.0*j30x+t11)*t4)*t5+(81.0*m10*t4+
. 729.0*t12)*t13+(-81.0*j30x*t4*t5+t14)*cos(q2)**2*cos(q3)**2
tt0=tt0+(-j10y*j30x*j30y+j10y*j30x*j30z+((t0+j30x*j30z)*m10+(-9.0*
. j10y*j30x+(-9.0*j10y+t1)*j30y+t10*j30z)*m30)*t3+((t1+t2+9.0*j30z)
. *m10*m30+(-81.0*j10y-81.0*j30x-81.0*j30y+81.0*j30z)*t4)*t5+(-81.0
. *m10*t4-729.0*t12)*t13+(-81.0*j30y*t4*t5+t14)*t15)*t9
inv(1,1)=-(t0+(t1+t2)*m30*t3-81.0*t4*t5+(j30x*j30y-j30x*j30z+(t6+
. t7-9.0*j30z)*m30*t3+81.0*t4*t5)*t9)//(tt0+((t11-81.0*j30z)*t4*t5+
. 729.0*t12*t13)*t15*t8**4)
t0=-j30x*j30y
t1=-9.0*j30x
t2=-9.0*j30y
t3=p**2
t4=m30**2
t5=p**4
t6=-81.0*t4*t5
t7=cos(q2)
t8=cos(q3)
t9=9.0*j30x
t10=9.0*j30y
t11=sin(q3)
t12=t11**2
t13=9.0*j10y+t9
t14=81.0*j30y
t15=m30**3
t16=p**6
t17=-729.0*t15*t16
t18=sin(q2)**2
tt0=j10y*j30x*j30y+(j30x*j30y*m10+(9.0*j10y*j30x+t13*j30y)*m30)*t3
. +((t9+t10)*m10*m30+(81.0*j10y+81.0*j30x+t14)*t4)*t5+(81.0*m10*t4+
. 729.0*t15)*t16+(-81.0*j30x*t4*t5+t17)*t7**2*t8**2
tt0=tt0+(-j10y*j30x*j30y+j10y*j30x*j30z+((t0+j30x*j30z)*m10+(-9.0*
. j10y*j30x+(-9.0*j10y+t1)*j30y+t13*j30z)*m30)*t3+((t1+t2+9.0*j30z)
. *m10*m30+(-81.0*j10y-81.0*j30x-81.0*j30y+81.0*j30z)*t4)*t5+(-81.0
. *m10*t4-729.0*t15)*t16+(-81.0*j30y*t4*t5+t17)*t18)*t12
inv(1,2)=(t0+(t1+t2)*m30*t3+t6+(-9.0*j30x*m30*t3+t6)*t7*t8+(j30x*
. j30y-j30x*j30z+(t9+t10-9.0*j30z)*m30*t3+81.0*t4*t5)*t12)//(tt0+((
. t14-81.0*j30z)*t4*t5+729.0*t15*t16)*t18*t11**4)
t0=p**2
t1=m30**2
t2=p**4
t3=sin(q2)
t4=sin(q3)
t5=9.0*j30y
t6=9.0*j30x
t7=9.0*j10y+t6
t8=81.0*j30y
t9=m30**3
t10=p**6
t11=-729.0*t9*t10
t12=-9.0*j30x
t13=t3**2
tt0=-j10y*j30x*j30y+j10y*j30x*j30z+((-j30x*j30y+j30x*j30z)*m10+(
. -9.0*j10y*j30x+(-9.0*j10y+t12)*j30y+t7*j30z)*m30)*t0+((t12-9.0*
. j30y+9.0*j30z)*m10*m30+(-81.0*j10y-81.0*j30x-81.0*j30y+81.0*j30z)
. *t1)*t2+(-81.0*m10*t1-729.0*t9)*t10
tt0=j10y*j30x*j30y+(j30x*j30y*m10+(9.0*j10y*j30x+t7*j30y)*m30)*t0+
. ((t6+t5)*m10*m30+(81.0*j10y+81.0*j30x+t8)*t1)*t2+(81.0*m10*t1+
. 729.0*t9)*t10+(-81.0*j30x*t1*t2+t11)*cos(q2)**2*cos(q3)**2+(tt0+(
. -81.0*j30y*t1*t2+t11)*t13)*t4**2
inv(1,3)=-((-9.0*j30y*m30*t0-81.0*t1*t2)*t3*t4+((t5-9.0*j30z)*m30*
. t0+81.0*t1*t2)*t3*t4**3)//(tt0+((t8-81.0*j30z)*t1*t2+729.0*t9*t10
. )*t13*t4**4)
t0=-j30x*j30y
t1=-9.0*j10y
t2=-9.0*j30y
t3=p**2
t4=m30**2
t5=p**4
t6=cos(q2)
t7=cos(q3)
t8=9.0*j30x
t9=9.0*j30y
t10=sin(q2)**2
t11=sin(q3)
t12=t11**2
t13=9.0*j10y+t8
t14=81.0*j30y
t15=m30**3
t16=p**6
t17=-729.0*t15*t16
t18=-9.0*j30x
tt1=-(-j10y*j30x+t0+(-j30x*m10+(t1-18.0*j30x+t2)*m30)*t3+(-9.0*m10
. *m30-162.0*t4)*t5+(-18.0*j30x*m30*t3-162.0*t4*t5)*t6*t7+(j30x*
. j30y-j30x*j30z+(t8+t9-9.0*j30z)*m30*t3+81.0*t4*t5+81.0*t4*t5*t10)
. *t12)
tt0=j10y*j30x*j30y+(j30x*j30y*m10+(9.0*j10y*j30x+t13*j30y)*m30)*t3
. +((t8+t9)*m10*m30+(81.0*j10y+81.0*j30x+t14)*t4)*t5+(81.0*m10*t4+
. 729.0*t15)*t16+(-81.0*j30x*t4*t5+t17)*t6**2*t7**2
tt0=tt0+(-j10y*j30x*j30y+j10y*j30x*j30z+((t0+j30x*j30z)*m10+(-9.0*
. j10y*j30x+(t1+t18)*j30y+t13*j30z)*m30)*t3+((t18+t2+9.0*j30z)*m10*
. m30+(-81.0*j10y-81.0*j30x-81.0*j30y+81.0*j30z)*t4)*t5+(-81.0*m10*
. t4-729.0*t15)*t16+(-81.0*j30y*t4*t5+t17)*t10)*t12
inv(2,2)=tt1//(tt0+((t14-81.0*j30z)*t4*t5+729.0*t15*t16)*t10*t11**
. 4)
t0=p**2
t1=m30**2
t2=p**4
t3=sin(q2)
t4=cos(q2)
t5=cos(q3)
t6=sin(q3)
t7=9.0*j30y
t8=9.0*j30x
t9=9.0*j10y+t8
t10=81.0*j30y
t11=m30**3
t12=p**6
t13=-729.0*t11*t12
t14=-9.0*j30x
t15=t3**2
tt0=-j10y*j30x*j30y+j10y*j30x*j30z+((-j30x*j30y+j30x*j30z)*m10+(
. -9.0*j10y*j30x+(-9.0*j10y+t14)*j30y+t9*j30z)*m30)*t0+((t14-9.0*
. j30y+9.0*j30z)*m10*m30+(-81.0*j10y-81.0*j30x-81.0*j30y+81.0*j30z)
. *t1)*t2+(-81.0*m10*t1-729.0*t11)*t12
tt0=j10y*j30x*j30y+(j30x*j30y*m10+(9.0*j10y*j30x+t9*j30y)*m30)*t0+
. ((t8+t7)*m10*m30+(81.0*j10y+81.0*j30x+t10)*t1)*t2+(81.0*m10*t1+
. 729.0*t11)*t12+(-81.0*j30x*t1*t2+t13)*t4**2*t5**2+(tt0+(-81.0*
. j30y*t1*t2+t13)*t15)*t6**2
inv(2,3)=(((-9.0*j30y*m30*t0-81.0*t1*t2)*t3-81.0*t1*t2*t4*t3*t5)*
. t6+((t7-9.0*j30z)*m30*t0+81.0*t1*t2)*t3*t6**3)//(tt0+((t10-81.0*
. j30z)*t1*t2+729.0*t11*t12)*t15*t6**4)
t0=-9.0*j10y
t1=-9.0*j30y
t2=p**2
t3=m30**2
t4=p**4
t5=cos(q2)**2
t6=cos(q3)**2
t7=9.0*j10y
t8=9.0*j30y
t9=sin(q3)
t10=t9**2
t11=9.0*j30x
t12=t7+t11
t13=81.0*j30y
t14=m30**3
t15=p**6
t16=-729.0*t14*t15
t17=-9.0*j30x
t18=sin(q2)**2
tt0=j10y*j30x*j30y+(j30x*j30y*m10+(9.0*j10y*j30x+t12*j30y)*m30)*t2
. +((t11+t8)*m10*m30+(81.0*j10y+81.0*j30x+t13)*t3)*t4+(81.0*m10*t3+
. 729.0*t14)*t15+(-81.0*j30x*t3*t4+t16)*t5*t6
tt0=tt0+(-j10y*j30x*j30y+j10y*j30x*j30z+((-j30x*j30y+j30x*j30z)*
. m10+(-9.0*j10y*j30x+(t0+t17)*j30y+t12*j30z)*m30)*t2+((t17+t1+9.0*
. j30z)*m10*m30+(-81.0*j10y-81.0*j30x-81.0*j30y+81.0*j30z)*t3)*t4+(
. -81.0*m10*t3-729.0*t14)*t15+(-81.0*j30y*t3*t4+t16)*t18)*t10
inv(3,3)=-(-j10y*j30y+(-j30y*m10+(t0+t1)*m30)*t2+(-9.0*m10*m30
. -81.0*t3)*t4+81.0*t3*t4*t5*t6+(j10y*j30y-j10y*j30z+((j30y-j30z)*
. m10+(t7+t8-9.0*j30z)*m30)*t2+(9.0*m10*m30+81.0*t3)*t4)*t10)//(tt0
. +((t13-81.0*j30z)*t3*t4+729.0*t14*t15)*t18*t9**4)
c
c --- Copy Entries Across Main Diagonals ---
c
do 25001 i=1,3
do 25002 j=i+1.0,3
mat(j,i)=mat(i,j)
inv(j,i)=inv(i,j)
25002 continue
25001 continue
--- NEW FILE: mcond.lisp ---
(setq m 10)
(gentran `( ((mdo) n ,m nil nil 100 nil
((mcond) ((mequal) k 3)
(($break))
t
$false)) ) nil)
--- NEW FILE: mdo.lisp ---
;; generating do loop from lisp
(setq n 10)
(gentran `(((mdo) i 1 nil nil ,n nil
((mdo) j ((mplus) i 1) nil nil ,n nil
((mprogn)
((msetq) ((x) j i) ((x) i j))
((msetq) ((y) j i) ((y) i j)))))) nil)
--- NEW FILE: rk.in ---
/* Runge-Kutta Method */
rungekutta(p1, p2, p, q, tt) :=
block( [k11, k12, k21, k22, k31, k32, k41, k42],
k11 : hh*p1,
k12 : hh*p2,
k21 : hh*ratsubst(q+k12/2, q,
ratsubst(p+k11/2, p,
ratsubst(tt+hh/2, tt, p1))),
k22 : hh*ratsubst(q+k12/2, q,
ratsubst(p+k11/2, p,
ratsubst(tt+hh/2, tt, p2))),
k31 : hh*ratsubst(q+k22/2, q,
ratsubst(p+k21/2, p,
ratsubst(tt+hh/2, tt, p1))),
k32 : hh*ratsubst(q+k22/2, q,
ratsubst(p+k21/2, p,
ratsubst(tt+hh/2, tt, p2))),
k41 : hh*ratsubst(q+k32, q,
ratsubst(p+k31, p,
ratsubst(tt+hh, tt, p1))),
k42 : hh*ratsubst(q+k32, q,
ratsubst(p+k31, p,
ratsubst(tt+hh, tt, p2))),
pn : ratsimp(p + (k11 + 2*k21 + 2*k31 + k41)/6),
qn : ratsimp(q + (k12 + 2*k22 + 2*k32 + k42)/6)
)$
--- NEW FILE: runge.mac ---
/* Enter parameters */
k : 1/(2*m)*p^2$
u : k0/2*q^2$
d : b/2*qdot$
h : k + u$
/* Perform Calculations */
batch("rk.in"); /* define Runge-Katta method */
batch("ham.in"); /* Hamiltonian Calculation */
gentranlang : fortran$
on(gentranfloat)$
/* Call Template Processor */
gentranin("runge.template", ["runge.f"]);
--- NEW FILE: runge.output ---
C fortran code generated from runge.template
program runge
implicit real (k,m)
c
c Input
c
write(6,*) 'Initial Value of p'
read(5,*) p
write(6,*) ' p = ', p
write(6,*) 'Initial Value of q'
read(5,*) q
write(6,*) ' q = ', q
write(6,*) 'Value of m'
read(5,*) m
write(6,*) ' m = ', m
write(6,*) 'Value of k0'
read(5,*) k0
write(6,*) ' k0 = ', k0
write(6,*) 'Value of b'
read(5,*) b
write(6,*) ' b = ', b
write(6,*) 'Step Size of t'
read(5,*) hh
write(6,*) ' Step Size of t = ', hh
write(6,*) 'Final Value of t'
read(5,*) tp
write(6,*) ' Final Value of t = ', tp
c
c Initialization
c
tt=0.0
write(9,*) ' h = p**2/2.0/m+k0*q**2/2.0'
write(9,*) ' d = b*qdot/2.0'
write(9,901) c
901 format(' c= ',e20.10)
write(9,910) tt, q, p
910 format(' ',3e20.10)
c
c Loop
c
25009 if (tt .ge. tf) goto 25010
pn=-1.0/2.0*b*hh-1.0/96.0*b*hh**5*k0**2/m**2+b*hh**3*k0/12.0/m
. +p-1.0/96.0*hh**6*k0**3*p/m**3+hh**4*k0**2*p/8.0/m**2-1.0/2.0
. *hh**2*k0*p/m-hh*k0*q-1.0/48.0*hh**5*k0**3*q/m**2+hh**3*k0**2
. *q/6.0/m
q=b*hh**8*k0**3/384.0/m**4-1.0/48.0*b*hh**6*k0**2/m**3+b*hh**4
. *k0/12.0/m**2-1.0/4.0*b*hh**2/m+hh**9*k0**4*p/384.0/m**5-1.0/
. 32.0*hh**7*k0**3*p/m**4+hh**5*k0**2*p/6.0/m**3-1.0/2.0*hh**3*
. k0*p/m**2+hh*p/m+q+hh**8*k0**4*q/192.0/m**4-1.0/24.0*hh**6*k0
. **3*q/m**3+hh**4*k0**2*q/6.0/m**2-1.0/2.0*hh**2*k0*q/m
p=pn
tt=tt+hh
write(9,910) tt, q, p
goto 25009
25010 continue
stop
end
--- NEW FILE: runge.template ---
C fortran code generated from runge.template
program runge
implicit real (k,m)
c
c Input
c
write(6,*) 'Initial Value of p'
read(5,*) p
write(6,*) ' p = ', p
write(6,*) 'Initial Value of q'
read(5,*) q
write(6,*) ' q = ', q
write(6,*) 'Value of m'
read(5,*) m
write(6,*) ' m = ', m
write(6,*) 'Value of k0'
read(5,*) k0
write(6,*) ' k0 = ', k0
write(6,*) 'Value of b'
read(5,*) b
write(6,*) ' b = ', b
write(6,*) 'Step Size of t'
read(5,*) hh
write(6,*) ' Step Size of t = ', hh
write(6,*) 'Final Value of t'
read(5,*) tp
write(6,*) ' Final Value of t = ', tp
c
c Initialization
c
tt=0.0
<<
gentran( literal(tab, "write(9,*) ' h = ", eval(h), "'", cr),
literal(tab, "write(9,*) ' d = ", eval(d), "'", cr) )$
>>
write(9,901) c
901 format(' c= ',e20.10)
write(9,910) tt, q, p
910 format(' ',3e20.10)
c
c Loop
c
<<
gentran( unless tt >= tf do
(
rsetq(pn, ev(pn,expand)),
rsetq(q, ev(qn,expand)),
p : pn,
tt : tt + hh,
literal(tab, "write(9,910) tt, q, p", cr)
) )$
>>
stop
end
--- NEW FILE: t.mac ---
gentran(x:y+t, if y=true then x:t else t:y);
--- NEW FILE: type.mac ---
gentranlang:fortran;
/* gentran knows about constants and implicit types in fortran */
gentran(i:5+j*2);
gentran(x:y/6-3);
/* types can also be explicitly declared */
gentran( type("implicit real*8", "a-h","o-z"),
type("real*8", m(4,4)),
for i:1 thru 4 do
for j:1 thru 4 do
if i=j
then m[i,j] : 1.0
else m[i,j] : 0.0,
type("integer", i,j));
gentran( type("double", x),
type("integer", y),
i:9*y,
x:7*x/u)$
--- NEW FILE: type.output ---
i=5+j*2
x=y/6.0-3.0
implicit real*8 (a-h,o-z)
real*8 m(4,4)
integer i,j
do 25009 i=1,4
do 25010 j=1,4
if ( .not. i .eq. j) goto 25011
m(i,j)=1.0
goto 25012
25011 continue
m(i,j)=0.0
25012 continue
25010 continue
25009 continue
double x
integer y
i=9*y
x=7.d0*x/u
|