(* This is the mathematica script that generates the
   text of ricci_calc.inc *)


<< /afs/ncsa.uiuc.edu/projects/genrel/home/mamiller/MATHEMATICA/Format.m;

dir1sub = {D[gxx[x,y,z],x] -> dexgxx,D[gxy[x,y,z],x] -> dexgxy,
           D[gxz[x,y,z],x] -> dexgxz,D[gyy[x,y,z],x] -> dexgyy,
           D[gyz[x,y,z],x] -> dexgyz,D[gzz[x,y,z],x] -> dexgzz,
           D[gxx[x,y,z],y] -> deygxx,D[gxy[x,y,z],y] -> deygxy,
           D[gxz[x,y,z],y] -> deygxz,D[gyy[x,y,z],y] -> deygyy,
           D[gyz[x,y,z],y] -> deygyz,D[gzz[x,y,z],y] -> deygzz,
           D[gxx[x,y,z],z] -> dezgxx,D[gxy[x,y,z],z] -> dezgxy,
           D[gxz[x,y,z],z] -> dezgxz,D[gyy[x,y,z],z] -> dezgyy,
           D[gyz[x,y,z],z] -> dezgyz,D[gzz[x,y,z],z] -> dezgzz};

dir2sub = {D[D[gxx[x,y,z],x],x] -> d2xxgxx,D[D[gxy[x,y,z],x],x] -> d2xxgxy,
           D[D[gxz[x,y,z],x],x] -> d2xxgxz,D[D[gyy[x,y,z],x],x] -> d2xxgyy,
           D[D[gyz[x,y,z],x],x] -> d2xxgyz,D[D[gzz[x,y,z],x],x] -> d2xxgzz,
           D[D[gxx[x,y,z],x],y] -> d2xygxx,D[D[gxy[x,y,z],x],y] -> d2xygxy,
           D[D[gxz[x,y,z],x],y] -> d2xygxz,D[D[gyy[x,y,z],x],y] -> d2xygyy,
           D[D[gyz[x,y,z],x],y] -> d2xygyz,D[D[gzz[x,y,z],x],y] -> d2xygzz,
           D[D[gxx[x,y,z],x],z] -> d2xzgxx,D[D[gxy[x,y,z],x],z] -> d2xzgxy,
           D[D[gxz[x,y,z],x],z] -> d2xzgxz,D[D[gyy[x,y,z],x],z] -> d2xzgyy,
           D[D[gyz[x,y,z],x],z] -> d2xzgyz,D[D[gzz[x,y,z],x],z] -> d2xzgzz,
           D[D[gxx[x,y,z],y],y] -> d2yygxx,D[D[gxy[x,y,z],y],y] -> d2yygxy,
           D[D[gxz[x,y,z],y],y] -> d2yygxz,D[D[gyy[x,y,z],y],y] -> d2yygyy,
           D[D[gyz[x,y,z],y],y] -> d2yygyz,D[D[gzz[x,y,z],y],y] -> d2yygzz,
           D[D[gxx[x,y,z],y],z] -> d2yzgxx,D[D[gxy[x,y,z],y],z] -> d2yzgxy,
           D[D[gxz[x,y,z],y],z] -> d2yzgxz,D[D[gyy[x,y,z],y],z] -> d2yzgyy,
           D[D[gyz[x,y,z],y],z] -> d2yzgyz,D[D[gzz[x,y,z],y],z] -> d2yzgzz,
           D[D[gxx[x,y,z],z],z] -> d2zzgxx,D[D[gxy[x,y,z],z],z] -> d2zzgxy,
           D[D[gxz[x,y,z],z],z] -> d2zzgxz,D[D[gyy[x,y,z],z],z] -> d2zzgyy,
           D[D[gyz[x,y,z],z],z] -> d2zzgyz,D[D[gzz[x,y,z],z],z] -> d2zzgzz};

dir3sub = {D[uxx[x,y,z],x] -> dxuxx,D[uxy[x,y,z],x] -> dxuxy,
           D[uxz[x,y,z],x] -> dxuxz,D[uyy[x,y,z],x] -> dxuyy,
           D[uyz[x,y,z],x] -> dxuyz,D[uzz[x,y,z],x] -> dxuzz,
           D[uxx[x,y,z],y] -> dyuxx,D[uxy[x,y,z],y] -> dyuxy,
           D[uxz[x,y,z],y] -> dyuxz,D[uyy[x,y,z],y] -> dyuyy,
           D[uyz[x,y,z],y] -> dyuyz,D[uzz[x,y,z],y] -> dyuzz,
           D[uxx[x,y,z],z] -> dzuxx,D[uxy[x,y,z],z] -> dzuxy,
           D[uxz[x,y,z],z] -> dzuxz,D[uyy[x,y,z],z] -> dzuyy,
           D[uyz[x,y,z],z] -> dzuyz,D[uzz[x,y,z],z] -> dzuzz};

outsub = {uxx[x,y,z] -> uxx,uxy[x,y,z] -> uxy,uxz[x,y,z] -> uxz,
          uyy[x,y,z] -> uyy,uyz[x,y,z] -> uyz,uzz[x,y,z] -> uzz};

coord = {x,y,z};
metric = {{gxx[x,y,z],gxy[x,y,z],gxz[x,y,z]},
          {gxy[x,y,z],gyy[x,y,z],gyz[x,y,z]},
          {gxz[x,y,z],gyz[x,y,z],gzz[x,y,z]}};
dmetric = Table[0,{i,3},{j,3},{k,3}];
For[i=1,i<=3,i++,
For[j=1,j<=3,j++,
For[k=1,k<=3,k++,
   dmetric[[i,j,k]] = D[metric[[j,k]],coord[[i]]];
];
];
];

invmetric = {{uxx[x,y,z],uxy[x,y,z],uxz[x,y,z]},
          {uxy[x,y,z],uyy[x,y,z],uyz[x,y,z]},
          {uxz[x,y,z],uyz[x,y,z],uzz[x,y,z]}};

Print["c     here are the assignments for dxuxx et al \n"];
Print[FortranAssign[dxuxx,(( - Sum[invmetric[[1,i]] invmetric[[1,j]] *
   dmetric[[1,i,j]],{i,3},{j,3}]) /. dir1sub) /. outsub]];
Print[FortranAssign[dxuxy,(( - Sum[invmetric[[1,i]] invmetric[[2,j]] *
   dmetric[[1,i,j]],{i,3},{j,3}]) /. dir1sub) /. outsub]];
Print[FortranAssign[dxuxz,(( - Sum[invmetric[[1,i]] invmetric[[3,j]] *
   dmetric[[1,i,j]],{i,3},{j,3}]) /. dir1sub) /. outsub]];
Print[FortranAssign[dxuyy,(( - Sum[invmetric[[2,i]] invmetric[[2,j]] *
   dmetric[[1,i,j]],{i,3},{j,3}]) /. dir1sub) /. outsub]];
Print[FortranAssign[dxuyz,(( - Sum[invmetric[[2,i]] invmetric[[3,j]] *
   dmetric[[1,i,j]],{i,3},{j,3}]) /. dir1sub) /. outsub]];
Print[FortranAssign[dxuzz,(( - Sum[invmetric[[3,i]] invmetric[[3,j]] *
   dmetric[[1,i,j]],{i,3},{j,3}]) /. dir1sub) /. outsub]];

Print[FortranAssign[dyuxx,(( - Sum[invmetric[[1,i]] invmetric[[1,j]] *
   dmetric[[2,i,j]],{i,3},{j,3}]) /. dir1sub) /. outsub]];
Print[FortranAssign[dyuxy,(( - Sum[invmetric[[1,i]] invmetric[[2,j]] *
   dmetric[[2,i,j]],{i,3},{j,3}]) /. dir1sub) /. outsub]];
Print[FortranAssign[dyuxz,(( - Sum[invmetric[[1,i]] invmetric[[3,j]] *
   dmetric[[2,i,j]],{i,3},{j,3}]) /. dir1sub) /. outsub]];
Print[FortranAssign[dyuyy,(( - Sum[invmetric[[2,i]] invmetric[[2,j]] *
   dmetric[[2,i,j]],{i,3},{j,3}]) /. dir1sub) /. outsub]];
Print[FortranAssign[dyuyz,(( - Sum[invmetric[[2,i]] invmetric[[3,j]] *
   dmetric[[2,i,j]],{i,3},{j,3}]) /. dir1sub) /. outsub]];
Print[FortranAssign[dyuzz,(( - Sum[invmetric[[3,i]] invmetric[[3,j]] *
   dmetric[[2,i,j]],{i,3},{j,3}]) /. dir1sub) /. outsub]];

Print[FortranAssign[dzuxx,(( - Sum[invmetric[[1,i]] invmetric[[1,j]] *
   dmetric[[3,i,j]],{i,3},{j,3}]) /. dir1sub) /. outsub]];
Print[FortranAssign[dzuxy,(( - Sum[invmetric[[1,i]] invmetric[[2,j]] *
   dmetric[[3,i,j]],{i,3},{j,3}]) /. dir1sub) /. outsub]];
Print[FortranAssign[dzuxz,(( - Sum[invmetric[[1,i]] invmetric[[3,j]] *
   dmetric[[3,i,j]],{i,3},{j,3}]) /. dir1sub) /. outsub]];
Print[FortranAssign[dzuyy,(( - Sum[invmetric[[2,i]] invmetric[[2,j]] *
   dmetric[[3,i,j]],{i,3},{j,3}]) /. dir1sub) /. outsub]];
Print[FortranAssign[dzuyz,(( - Sum[invmetric[[2,i]] invmetric[[3,j]] *
   dmetric[[3,i,j]],{i,3},{j,3}]) /. dir1sub) /. outsub]];
Print[FortranAssign[dzuzz,(( - Sum[invmetric[[3,i]] invmetric[[3,j]] *
   dmetric[[3,i,j]],{i,3},{j,3}]) /. dir1sub) /. outsub]];

chris[a_,b_,c_] := (1/2) Sum[invmetric[[i,c]] (dmetric[[a,b,i]] +
   dmetric[[b,a,i]] - dmetric[[i,a,b]]),{i,3}];

ricci[a_,b_] := Sum[ D[chris[a,b,i],coord[[i]]] - D[chris[i,b,i],coord[[a]]] +
   Sum[(chris[a,b,j]) (chris[j,i,i]) - 
       (chris[i,b,j]) (chris[j,a,i]),{j,3}],{i,3}];

Print["   \n"];
Print["c     here are the assignments for the x-x component of the 3-ricci\n"];
Print[FortranAssign[riccixx , 
   Simplify[((((ricci[1,1] /. dir1sub) /. dir2sub) /. dir3sub) /. outsub)]]];

Print["   \n"];
Print["c     here are the assignments for the x-y component of the 3-ricci\n"];
Print[FortranAssign[riccixy , 
   Simplify[((((ricci[1,2] /. dir1sub) /. dir2sub) /. dir3sub) /. outsub)]]];

Print["   \n"];
Print["c     here are the assignments for the x-z component of the 3-ricci\n"];
Print[FortranAssign[riccixz ,
   Simplify[((((ricci[1,3] /. dir1sub) /. dir2sub) /. dir3sub) /. outsub)]]];

Print["   \n"];
Print["c     here are the assignments for the y-y component of the 3-ricci\n"];
Print[FortranAssign[ricciyy ,
   Simplify[((((ricci[2,2] /. dir1sub) /. dir2sub) /. dir3sub) /. outsub)]]];

Print["   \n"];
Print["c     here are the assignments for the y-z component of the 3-ricci\n"];
Print[FortranAssign[ricciyz ,
   Simplify[((((ricci[2,3] /. dir1sub) /. dir2sub) /. dir3sub) /. outsub)]]];

Print["   \n"];
Print["c     here are the assignments for the z-z component of the 3-ricci\n"];
Print[FortranAssign[riccizz ,
   Simplify[((((ricci[3,3] /. dir1sub) /. dir2sub) /. dir3sub) /. outsub)]]];