/* Q with easy change at run time of Lpx and x[2], x[3] (QX) and homogenized W2
tetrad - W3 - (for spherical symmetric Lorentz transformations)*/
/*Includes also transformations LT12W, LT23W and LT13W between the Wx and
Rotations and Boosts around W3*/
/*Some out of service eliminated*/
#include <stdio.h>
#include <math.h>
#define DEF1
#define mrs (1 - rs/x[1])
#define mrs2 sqrt(mrs)
#define trab (double x[]){double z=0;
#define sa sin(x[2])
#define sb sin(x[3])
#define ca cos(x[2])
#define cb cos(x[3])
#define V(i,j) (*V[i][j])(double x[])
#define t x[0]
#define r x[1]
#define a x[2]
#define b x[3]
#define fjk for(j=0;j<=3;j++){for(k=0;k<=3;k++){
#define d double
#define rz return z;}
#define fz (double x[]){double z=0;int i;for(i=0;i<=3;i++){ z+=
#define PI 3.1415927
#define A (Lp1*(x[1]-0.03)) /*Lorentz parameters Boosts - r = 0.3*/
#define B2 (Lp2*(x[1]-0.03))
#define B3 (Lp3*(x[1]-0.03))
#define B (Lp4*(x[1]-0.03)) /*Rotations*/
#define R2 (Lp5*(x[1]-0.03))
#define R3 (Lp6*(x[1]-0.03))
/*#define A (Lp1*(x[1]-0.03)*(x[1]-0.03)) /*Lorentz parameters Boosts - r = 0.3*/
#define B2 (Lp2*(x[1]-0.03)*(x[1]-0.03))
#define B3 (Lp3*(x[1]-0.03)*(x[1]-0.03))
#define B (Lp4*(x[1]-0.03)*(x[1]-0.03)) /*Rotations*/
#define R2 (Lp5*(x[1]-0.03)*(x[1]-0.03))
#define R3 (Lp6*(x[1]-0.03)*(x[1]-0.03))*/
#define K (8*PI*7.425*1e-29) /*GN/c2 Einstein Gravitational Constant in cm/gr*/
#define Gn 7.425e-29
#define M (rs/(2*Gn))
#define En (-0.5*Gn*M*M/L1) /*Newton Gravitational Energy*/

/*Used symbols t,r,a,b,rs,mrs,mrs2,trab,sa,sb,ca,cb,x[],V[][],V(i,j),g00*/
/*g11,g22,g33,g[],W00,..,W[][],V00,..,VB00,..,VB200,..,V[][],VB[][],VB2[][]*/
/*VRab00,..,VRab[][],XYtrabLB00,..,LB200,..,LB[][],LB2[][],LRab00,..,LRab[][]*/
/*Lor,LorB1R1,fjk,Lp1,Lp2,ggNN,gg[],d,rz,fz,PI,A,B,K,GN,gr*/

#ifndef integral
d EG(d rs, d Lp1) { d z; z = 4*PI*rs/K + (8*PI*Lp1/K)*(sqrt(1 - rs/Lp1) - 1); return z;}
#endif
int J = 0; /*(0, , ,6) Index for LorX[J]()*/
d Lp1 = 0; /* Lorentz parameter - Boosts*/
d Lp2 = 0;
d Lp3 = 0;
d Lp4 = 0;/*Rotations*/
d Lp5 = 0;
d Lp6 = 0;
d *LP[7] = {&Lp1,&Lp1,&Lp2,&Lp3,&Lp4,&Lp5,&Lp6};


int ETA[4] = {1,-1,-1,-1};
double x[4]={7.00, 1.00, 1.00 , 1.00};/*double x[4]={7.00, 1.00, 1.5708 , 0.00};*/
double dx[4]={0.1, 0.0001, 0.0001, 0.0001};
double rs=1e-2;

/*Covariant metric*/
double g00 trab z= mrs ; return z;}
double g11 trab z= -(1/mrs) ; return z;}
double g22 trab z= - (r * r); return z;}
double g33 trab z= -(r * r * sa * sa); return z;}
double (*g[4])() = { g00, g11, g22, g33};
d dg(d x[]){d z; z = g00(x) * g11(x) * g22(x) * g33(x); rz

/*Contravariant metric*/
d gg00 trab z = 1/mrs; rz
d gg11 trab z = - mrs; rz
d gg22 trab z = - 1/(r *r); rz
d gg33 trab z = - 1/(r * r * sa * sa); rz
d (*gg[4])() = {gg00,gg11,gg22,gg33};

double F0(double x[]){return 0.00;}
double F1(d x[]){return 1.00;}

/*Christoffel Symbols*/
d G001 trab z = rs/(2 * r * (r - rs)); return z;}
d G010 trab z = G001(x); rz
d G100 trab z = G001(x) * mrs * mrs; return z;}
d G111 trab z = -G001(x); return z;}
d G122 trab z = -r * mrs; return z;}
d G133 trab z = -r * mrs * sa * sa; rz
d G212 trab z = 1/r; rz
d G221 trab z = G212(x); rz
d G233 trab z = -sa * ca; rz
d G313 trab z = 1/r; rz
d G331 trab z = G313(x); rz
d G323 trab z = ca / sa; rz
d G332 trab z = G323(x); rz
d (*G[4][4][4])() = {F0,G001,F0,F0,G010,F0,F0,F0,F0,F0,F0,F0,F0,F0,F0,F0,G100,F0,F0,F0,F0,G111,
F0,F0,F0,F0,G122,F0,F0,F0,F0,G133,F0,F0,F0,F0,F0,F0,G212,F0,F0,G221,F0,F0,F0,
F0,F0,G233,F0,F0,F0,F0,F0,F0,F0,G313,F0,F0,F0,G323,F0,G331,G332,F0};

#ifdef integral
d integral (d (*f)(), d idx, d L1, d L2){d z=0; d y[4]; d dy; int i; int n;
y[1] = L1;
dy = 0;
while(y[1] < L2){
y[1] += dy ;
dy = y[1]*y[1]*idx/2;
y[1] += dy ;
z+= f(y) * 2*dy;} rz
#endif

/* W or W1 is the solution tetrad*/

double W00 trab
z = 1/sqrt(mrs);
return z;}
double W11 trab z= mrs2 * sa * cb; return z;}
double W12 trab z= 1/r * ca * cb; return z;}
double W13 trab z= -(1/(r * sa) * sb); return z;}
double W21 trab z= mrs2 * sa * sb; return z;}
double W22 trab z= 1/r *ca * sb; return z;}
double W23 trab z= 1/(r * sa) * cb; return z;}
double W31 trab z= mrs2 * ca; return z;}
double W32 trab z= -(1/r * sa); return z;}
double (*W[4][4])() = {W00,F0,F0,F0,F0,W11,W12,W13,F0,W21,W22,W23,F0,W31,W32,F0};

d AnaliticGG00(void){d z; z = (1 - mrs2) * (1 - mrs2)/(r * r); rz
d AnaliticGG11(void){d z; z = (1 - mrs2) * (1 - mrs2)/(r * r) + (1 - 1/mrs2) * rs/(r*r*r); rz
d AnaliticGG22(void){d z; z = 0.5 * (1 - 1/mrs2) * rs/(r*r*r); rz
d (*AnaliticGG[4][4])(void)={AnaliticGG00,F0,F0,F0,F0,AnaliticGG11,F0,F0,F0,F0,AnaliticGG22,F0,F0,F0,F0,AnaliticGG22};

/*Contravariant Gravitational Energy Tensor in gr/cm3*/
/*d integral (d (*f)(), d idx, d L0, d L1)*/
d TG00(void){d z; z = - gg00(x) * AnaliticGG[0][0]() / K ; rz

/*New part*/
/*W2 is the tetrad oriented as the spherical coordinates*/
double W200 trab z= 1/mrs2 ; return z;}
double W211 trab z= mrs2 ; return z;}
double W223 trab z= 1/(r * sa) ; return z;}
double W232 trab z= -1/r ; return z;}
double (*W2[4][4])() = {W200,F0,F0,F0,F0,W211,F0,F0,F0,F0,F0,W223,F0,F0,W232,F0};
/*LT23W Lorentz Transformation from W2 to W3 - W3 is a homogenized W2 for spherical
symmetric rotations*/
/*d LT23W22 trab z = cos (b); return z;}
d LT23W33 trab z = cos (b); return z;}
d LT23W23 trab z = sin (b); return z;}
d LT23W32 trab z = -sin (b); return z;}
d (*LT23W[4][4])() = {F1,F0,F0,F0,F0,F1,F0,F0,F0,F0,LT23W22,LT23W23,F0,F0,LT23W32,LT23W33};*/
/*LT23WI inverse LT23W*/
/*d (*LT23WI[4][4])() = {F1,F0,F0,F0,F0,F1,F0,F0,F0,F0,LT23W22,LT23W32,F0,F0,LT23W23,LT23W33};*/
/*LT21W Lorentz transformation from W2 to W1*/
/*d LT21W11 trab z = sa*cb; return z;}
d LT21W12 trab z = -sb; return z;}
d LT21W13 trab z = -ca*cb; return z;}
d LT21W21 trab z = sa*sb; return z;}
d LT21W22 trab z = cb; return z;}
d LT21W23 trab z = -ca*sb; return z;}
d LT21W31 trab z = ca; return z;}
d LT21W33 trab z = sa; return z;}
d (*LT21W[4][4])() = {F1,F0,F0,F0,F0,LT21W11,LT21W12,LT21W13,F0,LT21W21,LT21W22,LT21W23,F0,LT21W31,F0,LT21W33};*/
/*LT21WI inverse LT21W*/
/*d (*LT21WI[4][4])() = {F1,F0,F0,F0,F0,LT21W11,LT21W21,LT21W31,F0,LT21W12,LT21W22,F0,F0,LT21W13,LT21W23,LT21W33};*/
/*LT31W Lorentz Transformation W3 to W1*/
d LT31W00 trab z = 1.00; return z;}
d LT31W11 trab z = sa*cb; return z;}
d LT31W12 trab z = -sb*cb*(1+ca); return z;}
d LT31W13 trab z = sb*sb-ca*cb*cb; return z;}
d LT31W21 trab z = sa*sb; return z;}
d LT31W22 trab z = cb*cb-ca*sb*sb; return z;}
d LT31W23 trab z = -sb*cb*(1+ca); return z;}
d LT31W31 trab z = ca; return z;}
d LT31W32 trab z = sa*sb; return z;}
d LT31W33 trab z = sa*cb; return z;}
d (*LT31W[4][4])() = {F1,F0,F0,F0,F0,LT31W11,LT31W12,LT31W13,F0,LT31W21,LT31W22,LT31W23,F0,LT31W31,LT31W32,LT31W33};
/*LT31WI inverse of LT31W*/
d (*LT31WI[4][4])() = {F1,F0,F0,F0,F0,LT31W11,LT31W21,LT31W31,F0,LT31W12,LT31W22,LT31W32,F0,LT31W13,LT31W23,LT31W33};
/*Boosts*/
double LB100 (double x[]){double z; z=cosh(A);return z;}
double LB101 (double x[]){double z; z=sinh(A); return z;}
double (*LB1[4][4])() = {LB100,LB101,F0,F0,LB101,LB100,F0,F0,F0,F0,F1,F0,F0,F0,F0,F1};

double LB200 (double x[]){double z; z=cosh(B2);return z;}
double LB202 (double x[]){double z; z=sinh(B2); return z;}
double (*LB2[4][4])() = {LB200,F0,LB202,F0,F0,F1,F0,F0,LB202,F0,LB200,F0,F0,F0,F0,F1};

double LB300 (double x[]){double z; z=cosh(B3);return z;}
double LB303 (double x[]){double z; z=sinh(B3); return z;}
double (*LB3[4][4])() = {LB300,F0,F0,LB303,F0,F1,F0,F0,F0,F0,F1,F0,LB303,F0,F0,LB300};

/*Rotations */
d LR122(d x[]){d z; z = cos(B);return z;}
d LR123(d x[]){d z; z = - sin(B);return z;}
d LR132(d x[]){d z; z = sin(B);return z;}
d (*LR1[4][4])() ={F1,F0,F0,F0,F0,F1,F0,F0,F0,F0,LR122,LR123,F0,F0,LR132,LR122};

d LR211(d x[]){d z; z = cos(R2);return z;}
d LR213(d x[]){d z; z = sin(R2);return z;}
d LR231(d x[]){d z; z = -sin(R2);return z;}
d (*LR2[4][4])() ={F1,F0,F0,F0,F0,LR211,F0,LR213,F0,F0,F1,F0,F0,LR231,F0,LR211};

d LR311(d x[]){d z; z = cos(R3);return z;}
d LR312(d x[]){d z; z = -sin(R3);return z;}
d LR321(d x[]){d z; z = sin(R3);return z;}
d (*LR3[4][4])() ={F1,F0,F0,F0,F0,LR311,LR312,F0,F0,LR321,LR311,F0,F0,F0,F0,F1};

d (*Lor)();
void NoLor(void){fputs("No Lorentz Transformation Defined\n",stdout);}
d LorB1(d x[],int i, int j){d z ; z = LB1[i][j](x); return z;}
d LorB2(d x[],int i, int j){d z ; z = LB2[i][j](x); return z;}
d LorB3(d x[],int i, int j){d z ; z = LB3[i][j](x); return z;}
d LorR1(d x[],int i, int j){d z ; z = LR1[i][j](x); return z;}
d LorR2(d x[],int i, int j){d z ; z = LR2[i][j](x); return z;}
d LorR3(d x[],int i, int j){d z ; z = LR3[i][j](x); return z;}
d LorB1R1(d x[], int i, int j){d z = 0; int k; for(k=0;k<=3;k++) z += LB1[i][k](x) * LR1[k][j](x);
return z;}
d LorB1B2(d x[], int i, int j){d z = 0; int k; for(k=0;k<=3;k++) z += LB1[i][k](x) * LB2[k][j](x);
return z;}
d LorR1R3(d x[], int i, int j){d z = 0; int k; for(k=0;k<=3;k++) z += LR1[i][k](x) * LR3[k][j](x); return z;}
d Lor0(d x[], int i, int j){d z = 0.0; if(i==j) z = 1.0; return z;};
d (*LorX[7])()={Lor0, LorB1, LorB2, LorB3, LorR1, LorR2, LorR3};
d LorLT31W(d x[],int i, int j){d z ; z = LT31W[i][j](x); return z;}
d LorLTWR (d x[],int i, int j){d z=0;int k, l;
for(k=0;k<=3;k++) for(l=0;l<=3;l++)
z += LT31W[i][k](x)*LorX[J](x,k,l)*LT31WI[l][j](x);
return z;}
/*d LorLT21W(d x[],int i, int j){d z ; z = LT21W[i][j](x); return z;}
d LorLT23W(d x[],int i, int j){d z ; z = LT23W[i][j](x); return z;}*/

/*V is the tetrad W after Lorentz transformation*/
d V00 fz Lor(x,0,i) * W[i][0](x);} rz
d V01 fz Lor(x,0,i) * W[i][1](x);} rz
d V02 fz Lor(x,0,i) * W[i][2](x);} rz
d V03 fz Lor(x,0,i) * W[i][3](x);} rz
d V10 fz Lor(x,1,i) * W[i][0](x);} rz
d V11 fz Lor(x,1,i) * W[i][1](x);} rz
d V12 fz Lor(x,1,i) * W[i][2](x);} rz
d V13 fz Lor(x,1,i) * W[i][3](x);} rz
d V20 fz Lor(x,2,i) * W[i][0](x);} rz
d V21 fz Lor(x,2,i) * W[i][1](x);} rz
d V22 fz Lor(x,2,i) * W[i][2](x);} rz
d V23 fz Lor(x,2,i) * W[i][3](x);} rz
d V30 fz Lor(x,3,i) * W[i][0](x);} rz
d V31 fz Lor(x,3,i) * W[i][1](x);} rz
d V32 fz Lor(x,3,i) * W[i][2](x);} rz
d V33 fz Lor(x,3,i) * W[i][3](x);} rz

d (*V[4][4])()={V00,V01,V02,V03,V10,V11,V12,V13,V20,V21,V22,V23,V30,V31,V32,V33};
/*End of New part*/
d iAnaliticGG00(d x[]){d z; z = (1 - mrs2) * (1 - mrs2)/(r * r); rz
/*d iTG00(d x[]){d z; z = - 4*PI*r*r*gg00(x) * iAnaliticGG00(x) / K ; rz*/
d iEG(d x[]){d z; z = - 4*PI* V00(x) * (1 - mrs2) * (1 - mrs2) / K ; rz
double VgV (double (*V[][4])(), double (*U[][4])(), double x[], int j, int k)
{
int i;
char buf[30];
double z = 0;
for (i=0; i<=3; i++)
z += V[j][i](x) * g[i](x) * U[k][i](x);
return z;
}

double deriv(double (*f)(double x[]), double x[], double dx[], int i)
{double z;
double y2[4];
double y1[4];
int m;
for (m=0;m<=3;m++) {y1[m]=x[m];y2[m]=x[m];}
y2[i] = y2[i] + (dx[i]/2);

y1[i] = y1[i] - (dx[i]/2);
z = (f(y2) - f(y1))/dx[i];
return z;}


double derivM(double (*contderV)(),double (*V[][4])(), double x[], double dx[],int c,int j,int k,int l)
{double z;
d p;
d q;
int i;
double y1[4];
double y2[4];
for(i=0;i<=3;i++) {y1[i]=x[i];y2[i]=x[i];}
y2[l]=y2[l]+(dx[l]/2);

y1[l]=y1[l]-(dx[l]/2);
z = ((p=contderV(V,y2,dx,c,j,k)) - (q=contderV(V,y1,dx,c,j,k)))/dx[l];
return z;}

double GM (double (*contderV)(),double (*V[][4])(), d x[], d dx[], int c,int j,int k,int l)
{d z=0;
int i;
for(i=0;i<=3;i++)
z += G[j][l][i](x) * contderV(V,x,dx,c,i,k);
return z;}

d MG (d (*contderV)(), d (*V[][4])(), d x[], d dx[], int c,int j,int k,int l)
{d z=0;
int i;
for(i=0;i<=3;i++)
z += G[i][k][l](x) * contderV(V,x,dx,c,j,i);
return z;}


double GV (double (*V[][4])(), double x[], int c, int i, int j)
{double z = 0;
int m;
for(m=0;m<=3;m++)
z += G[i][j][m](x) * V[c][m](x); return z;}

double contderV (double (*V[][4])(double x[]), double x[], double dx[], int c, int j, int k)
{double z;
z = deriv(V[c][j], x, dx, k) + GV (V, x, c, j, k); return z;}

double mixderMV(d (*contderV)(d (*V[][4])(d x[]),d x[],d dx[],int c,int j,int k),d (*V[][4])(), d x[],d dx[],int c,int j,int k,int l)
{ d z;
z =derivM (contderV,V,x,dx,c,j,k,l) + GM(contderV,V,x,dx,c,j,k,l) - MG(contderV,V,x,dx,c,j,k,l);
return z;}
/*Ricci tensor when RM = - RG*/
double Rjl (d (*mixderMV)(),d (*V[][4])(),d x[],d dx[],int j,int l)
{d z=0;int c;int k;
for(c=0;c<=3;c++){for(k=0;k<=3;k++){
z += V[c][k](x) * ETA[c] * (mixderMV(contderV,V,x,dx,c,j,k,l) - mixderMV(contderV,V,x,dx,c,j,l,k));}}
rz

/*Pijk =Vai ETAab Vbj;k (i and j contravariant)*/
d Pijk (d (*V[][4])(), d x[], d dx[], int i, int j, int k)
{d z = 0; int c;
for(c=0;c<=3;c++)z += V[c][i](x) * ETA[c] * contderV(V, x, dx, c, j, k);
rz
/*Pijkl = Pijk;l*/
d Pijkl (d (*V[][4])(), d x[], d dx[], int i, int j, int k, int l){d z=0;d y=0;int m;
for(m=0;m<=3;m++) {y += G[i][m][l](x) * Pijk(V, x, dx, m, j, k)
+ G[j][m][l](x) * Pijk(V,x,dx,i,m,k) - G[m][k][l](x) * Pijk(V,x,dx,i,j,m);}
z =derivM (Pijk,V,x,dx,i,j,k,l) + y; rz

/*GGij = Gravitational part of Einstein tensor*/
/*RGij = Pkik;j - Pkij;k*/
void GGij(d GG[][4], d (*V[][4])(), d x[], d dx[])
{d z = 0; d y = 0; int j; int k; int m;
fjk for(m=0;m<=3;m++) z+=Pijkl(V,x,dx,m,j,m,k) - Pijkl(V,x,dx,m,j,k,m);
GG[j][k]= z; z=0;}}
for(j=0;j<=3;j++) y += GG[j][j]/2;
for(j=0;j<=3;j++) GG[j][j] = GG[j][j] - y;
}
#ifndef FOUR_COMPONENTS
double TGG(d GG[][4], d x[], d dx[])
{ d y = 0; int j; for(j=0;j<=3;j++) y += GG[j][j]; return y;}
#endif
/* double contderV (double (*V[][4])(double x[]), double x[], double dx[], int c, int j, int k)*/
/* RM2lm = Pijm Pjli - Piji Pjlm = Vai;i ETAaa Val;m - Vai;m ETAaa Val;i */
double RM2lm(double (*V[][4])(), double x[], double dx[],int l, int m)
{double z=0;int j; int k;
fjk z+= contderV(V,x,dx,j,k,k) * ETA[j] * contderV(V,x,dx,j,l,m) -
contderV(V,x,dx,j,k,m) *ETA[j] * contderV(V,x,dx,j,l,k);}}return z;}

#ifdef DEF1
double *GMAT(double dLp1, double GM[][4])/*Non used RM2lm - 1/2 delta RM2*/
{
int y = 0; int z;int i;int j;int k;
for (j=0;j<=3;j++){
for (k=0;k<=3;k++){
GM[j][k] = RM2lm(V,x,dx,j,k);
if(j == k) y += GM[j][k];}}
for(i=0;i<=3;i++) GM[i][i] = GM[i][i] - y/2;return GM;
}

/*VgVj = Vck gkk Vek;j*/
double VgVj (double (*V[][4])(), d x[], d dx[], int c, int e,int j)
{d z = 0;int i;
for(i=0;i<=3;i++)
z += V[c][i](x) * g[i](x) * contderV(V,x,dx,e,i,j);
return z;}

/*SRdgVVjj = (sqrt(-|g|) Vck Vek;j),j)*/
double SRdgVVjj (d (*f)(), d x[], d dx[], int c, int e)
{ d z = 0; d X1[4]; d X2[4]; int j;
for(j=0;j<=3;j++)
X1[j]=X2[j]=x[j];
for(j=0;j<=3;j++)
{
X1[j] = X1[j] - dx[j]/2;
X2[j] = X2[j] + dx[j]/2;
z += (sqrt(- dg(X2)) * f(V,X2,dx,c,e,j) * gg[j](X2) - sqrt(- dg(X1)) * f(V,X1,dx,c,e,j) * gg[j](X1))/dx[j] ;
}
return z;}

/* VgVjgj = (Vck Vek;j);j*/
double VgVjgj (d (*f)(), d x[], d dx[],int c, int e)
{d y = 0; d z = 0; d X1[4]; d X2[4]; int j; int k;
for(j=0;j<=3;j++)
X1[j]=X2[j]=x[j];
for(j=0;j<=3;j++)
{
X1[j] = X1[j] - dx[j]/2;
X2[j] = X2[j] + dx[j]/2;
y += (f(V,X2,dx,c,e,j) - f(V,X1,dx,c,e,j))/dx[j] * gg[j](x) ;
for(k=0;k<=3;k++)
z +=G[k][j][j](x) * f(V,x,dx,c,e,k) * gg[j](x);
}
z = y - z; return z;}

#endif DEF1

/*double mixderMV(d (*contderV)(d (*V[][4])(d x[]),d x[],d dx[],int c,int j,int k),d (*V[][4])(), d x[],d dx[],int c,int j,int k,int l)*/

void main (void)
{
int c;
int i;
int j;
int k;
int m;
int n;
double z;
double y;
double GM[4][4];
double G1[4][4]; /*FOUR COMPONENTS of the EINSTEIN Tensor*/
double G2[4][4];
double G3[4][4];
double G4[4][4];
double TG1; /* Auxiliar for Gx[4][4]*/
double TG2;
double TG3;
double TG4;
char ch[4];
static char buf[30];
char *buf2 = "?";
char *buf3 = "W1 enabled";
d GG[4][4];
d L1 = 7e10;
d L2 = 1e12;
d idx = 1e-15;
int graph = 0;/*Flag */
int LT = 1; /*Boost or Rotation (1,,,6) for graph*/


/*Save original tetrads*/
d (*W1[4][4])() ;
d (*Vx[4][4])() ;
fjk W1[j][k] = W[j][k];}}
fjk Vx[j][k] = V[j][k];}}
a1:
/*Ask for parameter values*/
graph = 0;
fputs("graph? (Y/N)",stdout);
gets(buf);
if (buf[0] == 'Y' || buf [0] == 'y') graph = 1;
fputs("rs = ",stdout);
fputs(gcvt(rs, 12, buf),stdout);fputs(" rs = ",stdout);
gets(buf);
if (buf[0] != '')
rs = atof(buf);

fputs("r = ",stdout);
fputs(gcvt(r, 12, buf),stdout);fputs(" r = ",stdout);
gets(buf);
if (buf[0] != '')
r = atof(buf);

fputs("dr = ",stdout);
fputs(gcvt(dx[1], 12, buf),stdout);fputs(" dr = ",stdout);
gets(buf);
if (buf[0] != '')
dx[1] = atof(buf);

fputs("dx[2] = ",stdout);fputs(gcvt(dx[2],12,buf),stdout);
fputs(" dx[2] = ",stdout);
gets(buf);
if (buf[0] != '')
dx[2] = atof(buf);

fputs("dx[3] = ",stdout);
fputs(gcvt(dx[3], 12, buf),stdout); fputs(" dx[3] = ",stdout);
gets(buf);
if (buf[0] != '')
dx[3] = atof(buf);

fputs("First Lorentz parameter = ", stdout); fputs(gcvt(Lp1, 12, buf),stdout);
fputs(" Lp1 = ",stdout);gets(buf); if(buf[0] != '') Lp1 = atof(buf);

fputs("Second Lorentz parameter = ", stdout); fputs(gcvt(Lp2, 12, buf),stdout);
fputs(" Lp2 = ",stdout);gets(buf); if(buf[0] != '') Lp2 = atof(buf);

fputs("Third Lorentz parameter = ", stdout); fputs(gcvt(Lp3, 12, buf),stdout);
fputs(" Lp3 = ",stdout);gets(buf); if(buf[0] != '') Lp3 = atof(buf);

fputs("Fourth Lorentz parameter = ", stdout); fputs(gcvt(Lp4, 12, buf),stdout);
fputs(" Lp4 = ",stdout);gets(buf); if(buf[0] != '') Lp4 = atof(buf);

fputs("Fifth Lorentz parameter = ", stdout); fputs(gcvt(Lp5, 12, buf),stdout);
fputs(" Lp5 = ",stdout);gets(buf); if(buf[0] != '') Lp5 = atof(buf);

fputs("6th Lorentz parameter = ", stdout); fputs(gcvt(Lp6, 12, buf),stdout);
fputs(" Lp6 = ",stdout);gets(buf); if(buf[0] != '') Lp6 = atof(buf);

fputs("Low limit of integration = ", stdout); fputs(gcvt(L1, 12, buf),stdout);
fputs(" L1 = ",stdout);gets(buf); if(buf[0] != '') L1 = atof(buf);

fputs("High limit of integration = ", stdout); fputs(gcvt(L2, 12, buf),stdout);
fputs(" L2 = ",stdout);gets(buf); if(buf[0] != '') L2 = atof(buf);

fputs("Integral increment = ", stdout); fputs(gcvt(idx, 12, buf),stdout);
fputs(" idx = ",stdout);gets(buf); if(buf[0] != '') idx = atof(buf);

/*Print parameter values*/
fputs("Lorentz transformation = ",stdout);fputs(buf2,stdout);fputs("\n",stdout);
fputs("rs = ",stdout); fputs(gcvt(rs,12,buf),stdout);
fputs(" r = ",stdout); fputs(gcvt(r,12,buf),stdout);
fputs(" dx[1] = ", stdout);fputs(gcvt(dx[1],12, buf), stdout);
fputs(" dx[2] = ", stdout);fputs(gcvt(dx[2],12, buf), stdout);
fputs(" dx[3] = ", stdout); fputs(gcvt(dx[3], 12, buf), stdout);fputs("\n",stdout);
fputs(" Lp1 = ",stdout); fputs(gcvt(Lp1,12,buf),stdout);
fputs(" Lp2 = ",stdout); fputs(gcvt(Lp2,12,buf),stdout);
fputs(" Lp3 = ",stdout); fputs(gcvt(Lp3,12,buf),stdout);
fputs(" Lp4 = ",stdout); fputs(gcvt(Lp4,12,buf),stdout);
fputs(" Lp5 = ",stdout); fputs(gcvt(Lp5,12,buf),stdout);
fputs(" Lp6 = ",stdout); fputs(gcvt(Lp6,12,buf),stdout);
fputs(" L1 = ",stdout); fputs(gcvt(L1,12,buf),stdout);
fputs(" L2 = ",stdout); fputs(gcvt(L2,12,buf),stdout);
fputs(" idx = ",stdout); fputs(gcvt(idx,12,buf),stdout);
fputs("\n",stdout);


fputs("rs = ",stdprn); fputs(gcvt(rs,12,buf),stdprn);
fputs(" r = ",stdprn); fputs(gcvt(r,12,buf),stdprn);
fputs(" dx[1] = ", stdprn);fputs(gcvt(dx[1],12, buf), stdprn);
fputs(" dx[2] = ", stdprn);fputs(gcvt(dx[2],12, buf), stdprn);
fputs(" dx[3] = ", stdprn); fputs(gcvt(dx[3], 12, buf), stdprn);
fputs("\n",stdprn);
fputs(" Lp1 = ",stdprn); fputs(gcvt(Lp1,12,buf),stdprn);
fputs(" Lp2 = ",stdprn); fputs(gcvt(Lp2,12,buf),stdprn);
fputs(" Lp3 = ",stdprn); fputs(gcvt(Lp3,12,buf),stdprn);
fputs(" Lp4 = ",stdprn); fputs(gcvt(Lp4,12,buf),stdprn);
fputs(" Lp5 = ",stdprn); fputs(gcvt(Lp5,12,buf),stdprn);
fputs(" Lp6 = ",stdprn); fputs(gcvt(Lp6,12,buf),stdprn);
fputs(" L1 = ",stdprn); fputs(gcvt(L1,12,buf),stdprn);
fputs(" L2 = ",stdprn); fputs(gcvt(L2,12,buf),stdprn);
fputs(" idx = ",stdprn); fputs(gcvt(idx,12,buf),stdprn);
fputs("\n",stdprn);

a0:
/*Menu of Lorentz transformations*/
fputs("Lorentz transformation = ",stdprn);fputs(buf2,stdprn);fputs("\n",stdprn);
fputs("Lorentz transformation = ",stdout);fputs(buf2,stdout);fputs("\n",stdout);
fputs("\nN No Lorentz 1 Boost1 2 Boost2 3 Boost3 4 Rot1 5 Rot2 6 Rot3 ENTER O.K.\n",stdout);
fputs("7 LT31W 8 LTWR 12 Boost1xBoost2 14 Rot1xRot3\n",stdout);
gets(buf);
if(buf[0] == '1'&& buf[1] == '4') goto R1R3;
if(buf[0] == '1'&& buf[1] == '2') goto BB;
if(buf[0] == '1') goto b1;
if(buf[0] == '2') goto b2;
if(buf[0] == '3') goto b3;
if(buf[0] == '4') goto r1;
if(buf[0] == '5') goto r2;
if(buf[0] == '6') goto r3;
if(buf[0] == '7') goto LTW31;
if(buf[0] == '8') goto LTWR;
if(buf[0] == 'N') goto nl;
if(buf[0] == '' && buf2[0] == '?') goto a0;
if(buf[0] == '') goto a00;
goto a0;
a00:
/*General Menu - The tetrads are contravariant and the tensors mixed*/
/*The Lorentz indeces are covariant*/
fputs("\n1 Paramters 2 V 3 VgV 4 Vcj,k 5 Vcj;k 6 Vcj;kl 7 Rjk 8 GMjk\n",stdout);
fputs("L Lorentz Matrix LT Lorentz Transformation W1 or W2 Tetrad 9 (VaiVbi;j);j\n",stdout);
fputs("10 GGjk 11 AnaliticGGjk 12 Contravariant TG00 13 Integral 14 FOUR COMPONENTS\n15 EXIT",stdout);
gets(buf);

if(buf[0] == 'L' && buf[1] == 'T') goto a0;
if(buf[0] == 'W' && buf[1] == '1') goto aW1;
if(buf[0] == 'W' && buf[1] == '2') goto aW2;
if(buf[0] == '1' && buf[1] == '0') goto a10;
if(buf[0] == '1' && buf[1] == '1') goto a11;
if(buf[0] == '1' && buf[1] == '2') goto a12;
if(buf[0] == '1' && buf[1] == '3') goto a13;
if(buf[0] == '1' && buf[1] == '4') goto a14;
if(buf[0] == '1' && buf[1] == '5') goto a15;
if(buf[0] == '1') goto a1;
if(buf[0] == '2') goto a2;
if(buf[0] == '3') goto a3;
if(buf[0] == '4') goto a4;
if(buf[0] == '5') goto a5;
if(buf[0] == '6') goto a6;
if(buf[0] == '7') goto a7;
if(buf[0] == '8') goto a8;
if(buf[0] == '9') goto a9;
if(buf[0] == 'L') goto L;
goto a00;
g0:
fputs("LT (0 - 6 Lorentz transformation = ",stdout);fputs(itoa(J,buf,10),stdout);
fputs(" LT = ",stdout);gets(buf);
if(buf[0] != '') LT = J = atoi(buf);
fputs("\n*Lorentz Paraneter[LT] = ",stdout);fputs(gcvt(*LP[LT],12,buf),stdout);
fputs(" *LP[LT] = ",stdout);gets(buf);
if(buf[0] != '') *LP[LT] = atof(buf);
fputs("\n*LP[",stdout);fputs(itoa(LT, buf, 10), stdout);fputs("] =",stdout);
fputs("\n*LP[",stdprn);fputs(itoa(LT, buf, 10), stdprn);fputs("] =",stdprn);
fputs(gcvt(*LP[LT],12,buf),stdout);fputs("\n",stdout);
fputs(buf,stdprn);fputs(" ",stdprn);
fputs("x[2] = ", stdout); fputs(gcvt(x[2],12,buf), stdout);
fputs(" x[2] = ", stdout); gets(buf);
if(buf[0] != '') x[2] = atof(buf);fputs("\n",stdout);
fputs("x[2] = ", stdprn); fputs(gcvt(x[2],12,buf),stdprn);fputs(" ",stdprn);
fputs("x[3] = ", stdout); fputs(gcvt(x[3],12,buf), stdout);
fputs(" x[3] = ", stdout); gets(buf);
if(buf[0] != '') x[3] = atof(buf);fputs("\n",stdout);
fputs("x[3] = ", stdprn); fputs(gcvt(x[3],12,buf),stdprn);fputs("\n",stdprn);
goto a00;

aW1:
fjk W[j][k] = W1[j][k];}}
buf3 = "W1 enabled";
buf2 = "?";
fputs(buf3,stdout);fputs("\n",stdout);
fputs(buf3,stdprn);fputs("\n",stdprn);
fputs("The tetrad is = \n",stdout);
fputs("The tetrad is = \n",stdprn);
fjk z = W[j][k](x);
fputs(gcvt(z, 12, buf), stdout);fputs(" ",stdout);
fputs(buf, stdprn);fputs(" ",stdprn);}
fputs("\n",stdout);fputs("\n",stdprn);}
fputs("Please indicate Lorentz Transformation\n",stdout);
goto a0;
aW2:
fjk W[j][k] = W2[j][k];}}
buf3 = "W2 enabled";
buf2 = "?";
fputs(buf3,stdout);fputs("\n",stdout);
fputs(buf3,stdprn);fputs("\n",stdprn);
fputs("The tetrad is = \n",stdout);
fputs("The tetrad is = \n",stdprn);
fjk z = W[j][k](x);
fputs(gcvt(z, 12, buf), stdout);fputs(" ",stdout);
fputs(buf, stdprn);fputs(" ",stdprn);}
fputs("\n",stdout);fputs("\n",stdprn);}
fputs("Please indicate Lorentz Transformation\n",stdout);
goto a0;
R1R3:
fjk V[j][k] = Vx[j][k];}}/*Return to V original value*/
Lor = LorR1R3;LT = 4;/*For LP[LT]*/
buf2 = "Rotation1 x Rotation3";
fputs(buf2,stdout);fputs("\n",stdout);
fputs(buf2,stdprn);fputs("\n",stdprn);
fputs("Lor = \n",stdout);
fputs("Lor = \n",stdprn);
fjk z = Lor(x,j,k);
fputs(gcvt(z, 12, buf), stdout);fputs(" ",stdout);
fputs(buf, stdprn);fputs(" ",stdprn);}
fputs("\n",stdout);fputs("\n",stdprn);}
goto a00;
LTWR:
buf2 = "LTWR\n"; LT = J;
fjk V[j][k] = Vx[j][k];}}/*Return to V original value*/
fputs(buf2,stdout);
fputs(buf2,stdprn);
Lor = LorLTWR;
goto a00;

LTW31:
buf2 = "LT31W\n";
fjk V[j][k] = Vx[j][k];}}/*Return to V original value*/
fputs(buf2,stdout);
fputs(buf2,stdprn);
Lor = LorLT31W;
goto a00;
b1:
buf2 = "Boost1";LT =1;
fjk V[j][k] = Vx[j][k];}}/*Return to V original value*/
fputs(buf2,stdout);
fputs(buf2,stdprn);
Lor = LorB1;
goto a00;
b2:
buf2 = "Boost2";LT = 2;
fjk V[j][k] = Vx[j][k];}}/*Return to V original value*/
fputs(buf2,stdout);fputs("\n",stdout);
fputs(buf2,stdprn);fputs("\n",stdprn);
Lor = LorB2;
goto a00;
b3:
buf2 = "Boost3";LT = 3;
fjk V[j][k] = Vx[j][k];}}/*Return to V original value*/
fputs(buf2,stdout);fputs("\n",stdout);
fputs(buf2,stdprn);fputs("\n",stdprn);
Lor = LorB3;
goto a00;
BB:
buf2 = "Boost1 x Boost2";
fjk V[j][k] = Vx[j][k];}}/*Return to V original value*/
fputs(buf2,stdout);fputs("\n",stdout);
fputs(buf2,stdprn);fputs("\n",stdprn);
Lor = LorB1B2;
goto a00;
r1:
buf2 = "Rotation on x1"; LT = 4;
fjk V[j][k] = Vx[j][k];}}/*Return to V original value*/
fputs(buf2,stdout);fputs("\n",stdout);
fputs(buf2,stdprn);fputs("\n",stdprn);
Lor = LorR1;
goto a00;
r2:
buf2 = "Rotation on x2";LT = 5;
fjk V[j][k] = Vx[j][k];}}/*Return to V original value*/
fputs(buf2,stdout);fputs("\n",stdout);
fputs(buf2,stdprn);fputs("\n",stdprn);
Lor = LorR2;
goto a00;
r3:
fjk V[j][k] = Vx[j][k];}}/*Return to V original value*/
buf2 = "Rotation on x3";LT = 6;
fputs(buf2,stdout);fputs("\n",stdout);
fputs(buf2,stdprn);fputs("\n",stdprn);
Lor = LorR3;
goto a00;
nl:
Lor = NoLor;
fjk V[j][k] = W[j][k];}}
buf2 = "Lorentz transformation disabled";
fputs(buf2,stdout);fputs("\n",stdout);
fputs(buf2,stdprn);fputs("\n",stdprn);
goto a00;

L:
if(Lor == NoLor){NoLor(); goto a00;}
fputs("L = \n",stdout);
fputs("L = \n",stdprn);
fjk z = Lor(x,j,k);
fputs(gcvt(z, 12, buf), stdout);fputs(" ",stdout);
fputs(buf, stdprn);fputs(" ",stdprn);}
fputs("\n",stdout);fputs("\n",stdprn);}
goto a00;

a2:
fputs(" k\n", stdout);
fputs("V = \n", stdout);
fputs(" a\n", stdout);
for (j=0;j<=3;j++){
for (k=0;k<=3;k++){
z = V[j][k](x);
fputs(gcvt(z, 12, buf), stdout);
fputs(" ", stdout);
}
fputs("\n", stdout);
}
goto a00;
a3:
#ifdef DEF1
fputs("\n VgV =\n", stdout);
for (j=0;j<=3;j++){
for (k=0;k<=3;k++){
z = VgV(V, V, x, j, k);

fputs(gcvt(z, 12, buf), stdout);
fputs(" ", stdout);
}
fputs("\n", stdout);
}
#endif DEF1
goto a00;
a4:
fputs("derivada Vcj,k = \n", stdout);

for(c=0;c<=3;c++)
{fputs("c = ", stdout); fputs(itoa(c, buf, 10),stdout);fputs("\n", stdout);
fjk z = deriv(V[c][j], x, dx, k);
fputs(gcvt(z, 12, buf), stdout); fputs(" ",stdout);}
fputs("\n",stdout);}}
goto a00;
a5:
for(c=0;c<=3;c++)
{
fputs(itoa(c, buf, 10), stdout);
fputs("covariant derivada Vcj;k = \n", stdout);
fjk z = contderV(V, x, dx, c, j, k);
fputs(gcvt(z, 12, buf), stdout); fputs(" ",stdout);}
fputs("\n",stdout);}
}
goto a00;
a6:
fputs("c = ",stdout); fputs(itoa (c, buf, 10), stdout);
gets(buf);
if (buf[0] != '')
c = atoi (buf);
fputs("c = ",stdout);
for(i=0;i<=3;i++)
{
fputs(itoa(c, buf, 10), stdout);fputs(", ",stdout);
fputs(itoa(i, buf, 10), stdout);
fputs(" Vcj;kl = \n", stdout);
fjk z = mixderMV(contderV, V, x, dx, c, i, j, k);
fputs(gcvt(z, 12, buf), stdout); fputs(" ",stdout);}
fputs("\n",stdout);}
}
goto a00;
/*Fuera de Servicio*/
#ifdef DEF2
fputs("\nGijk = \n ",stdout);
for(i=0;i<=3;i++){fjk if((z=G[i][j][k](x)) != 0) {c= i*100+j*10+k;fputs(itoa(c,buf,10),stdout);
fputs(" Gijk = ",stdout);fputs(gcvt(G[i][j][k](x),12,buf),stdout);
fputs(" ",stdout);}}}
fputs("\n",stdout);}

goto a0;
#endif DEF2
/*Fin de Fuera de Servicio*/
a7:
fputs(" j \n", stdout);
fputs("R = \n", stdout);
fputs(" l \n", stdout);
for (j=0;j<=3;j++){
for (k=0;k<=3;k++){
z = Rjl(mixderMV,V,x,dx,j,k);
fputs(gcvt(z, 12, buf), stdout);
fputs(" ", stdout);
}
fputs("\n", stdout);
}
goto a00;
a8:
y = 0;
for (j=0;j<=3;j++){
for (k=0;k<=3;k++){
z = RM2lm(V,x,dx,j,k);

GM[j][k] = z;
if(j == k) y += z;
}}
for(i=0;i<=3;i++) GM[i][i] = GM[i][i] - y/2;

fputs(" j \n",stdout);
fputs("GM =\n",stdout);
fputs(" k\n",stdout);
fputs("GM =\n",stdprn);
fjk fputs(gcvt(GM[j][k],12,buf),stdout); fputs(" ",stdout);
fputs(buf,stdprn);fputs(" ",stdprn);}
fputs("\n",stdout);
fputs("\n",stdprn);}
if(graph == 1) goto g0;
goto a00;
a9:
/*double VgVj (double (*V[][4])(), d x[], d dx[], int c, int e,int j)*/
/* double VgVjgj (d (*f)(), d x[], d dx[],int c, int e)*/
/*double SRdgVVjj (d (*V[][4])(), d x[], d dx[], int c, int e)*/
#ifdef DEF1
fputs("(VajVbj;k);k = \n",stdout);
fputs("(VajVbj;k);k = \n",stdprn);
fjk fputs(gcvt(VgVjgj(VgVj,x,dx,j,k),12,buf), stdout);fputs(" ",stdout);
fputs(buf, stdprn);fputs(" ",stdprn);}
fputs("\n",stdout);
fputs("\n",stdprn);}

fputs("(sqrt(-g)VajVbj;k),k = \n",stdout);
fputs("(sqrt(-g)VajVbj;k),k = \n",stdprn);
fjk fputs(gcvt(SRdgVVjj(VgVj,x,dx,j,k),12,buf), stdout);fputs(" ",stdout);
fputs(buf, stdprn);fputs(" ",stdprn);}
fputs("\n",stdout);
fputs("\n",stdprn);}
#endif DEF1
goto a00;
a10:
/*GGij = Gravitational part of Einstein Tensor*/
/*RGij = Pkjk;l - Pkjl;k*/
/*void GGij(d GG[][4], d (*V[][4])(), d x[], d dx[])*/
GGij(GG,V,x,dx);
fputs(" i\n",stdout);
fputs("GG = \n",stdout);
fputs(" j \n",stdout);
fputs("GGij = \n",stdprn);
fjk fputs(gcvt(GG[j][k],12,buf), stdout);fputs(" ",stdout);
fputs(buf, stdprn);fputs(" ",stdprn);}
fputs("\n",stdout);
fputs("\n",stdprn);}
#ifndef FOUR_COMPONENTS
fputs("Trace of GG = TGG = ", stdout);
fputs("Trace of GG = TGG = ", stdprn);
fputs(gcvt(TGG(GG, x, dx),12,buf), stdout);
fputs(buf, stdprn);
fputs("\n",stdout);
fputs("\n",stdprn);
#endif
if(graph == 1) goto g0;
goto a00;
a11:
fputs(" i\n",stdout);
fputs("AnaliticGG =\n",stdout);
fputs(" j\n",stdout);
fputs("AnaliticGGij = \n",stdprn);
fjk fputs(gcvt(AnaliticGG[j][k](),12,buf), stdout);fputs(" ",stdout);
fputs(buf, stdprn);fputs(" ",stdprn);}
fputs("\n",stdout);
fputs("\n",stdprn);}
goto a00;
a12:
fputs("If r and rs are given in cm. TG00 results in gr/cm3\n",stdout);

fputs("K = ",stdout);
fputs("K = ",stdprn);
fputs(gcvt(K,12,buf), stdout);fputs(" ",stdout);
fputs(buf, stdprn);fputs(" ",stdprn);
fputs("\n",stdout);
fputs("\n",stdprn);

fputs("Contravariant GG00 = ",stdout);
fputs("Contravariant GG00 = ",stdprn);
fputs(gcvt(gg00(x)*AnaliticGG00(),12,buf), stdout);fputs(" ",stdout);
fputs(buf, stdprn);fputs(" ",stdprn);
fputs("\n",stdout);
fputs("\n",stdprn);

fputs("Contravariant TG00 = ",stdout);
fputs("Contravariant TG00 = ",stdprn);
fputs(gcvt(TG00(),12,buf), stdout);fputs(" ",stdout);
fputs(buf, stdprn);fputs(" ",stdprn);
fputs("\n",stdout);
fputs("\n",stdprn);
goto a00;
a13:
/*Integral*/
/*d iTG00(d x[]){d z; z = - 4*PI*r*r*gg00(x) * iAnaliticGG00(x) / K ; rz*/
/*d iAnaliticGG00(d x[]){d z; z = (1 - mrs2) * (1 - mrs2)/(r * r); rz*/
/*d integral (d (*f)(), d idx, d L1, d L2){d z=0; d y[4]; d dy; int i; int n;*/
/* y[1] = L1;*/
/* while(y[1] < L2){*/
/* y[1] += dy = L1*L1*idx/2;*/
/* z+= f(y) * 2*dy;} rz*/

fputs("Integral of EG = sqrt(-g) V00 TG00 = ",stdout);
fputs("Integral of EG = sqrt(-g) V00 TG00 = ",stdprn);
#ifndef integral
fputs("Analitical Calculation = 4PIrs/k + 8PIL1/k (sqrt(1 - rs/L1) - 1) =\n",stdout);
#endif
#ifdef integral
fputs(gcvt(integral(iEG(),idx,L1,L2),12,buf), stdout);fputs(" ",stdout);
#endif
#ifndef integral
fputs(gcvt(EG(rs,L1),12,buf), stdout);fputs(" ",stdout);
#endif
fputs(buf, stdprn);fputs(" ",stdprn);
fputs("\n",stdout);
fputs("\n",stdprn);
fputs("E/En = ",stdout);fputs("E/En = ",stdprn);
fputs(gcvt(EG(rs,L1)/En,12,buf), stdout);
fputs(gcvt(EG(rs,L1)/En,12,buf), stdprn);
fputs("\n",stdout);
fputs("\n",stdprn);
goto a00;
a14:
/*Pijk =Vai ETAab Vbj;k (i and j contravariant)*/
/*d Pijk (d (*V[][4])(), d x[], d dx[], int i, int j, int k)*/
/*Pijkl = Pijk;l*/
/*d Pijkl (d (*V[][4])(), d x[], d dx[], int i, int j, int k, int l)*/
#ifdef FOUR_COMPONENTS
/* FOUR_COMPONENTS of EINSTEIN TENSOR - m contravariant, n covariant */
fputs("FOUR COMPONENTS of the EINSTEIN TENSOR",stdout);fputs("\n",stdout);
fputs("FOUR COMPONENTS of the EINSTEIN TENSOR",stdprn);fputs("\n",stdprn);

for (m=0;m<=3;m++){
for (n=0;n<=3;n++){
G1[m][n] = G2[m][n] = G3[m][n] = G4[m][n] = 0;}}
TG1 = TG2 = TG3 = TG4 = 0;

for (m=0;m<=3;m++){
for (n=0;n<=3;n++){
for (j=0;j<=3;j++){
G1[m][n] += Pijkl(V, x, dx, j, m, j, n);}}}

for (j=0;j<=3;j++){
TG1 += G1[j][j];}

for (j=0;j<=3;j++){
G1[j][j] = G1[j][j] - TG1/2;}

for (m=0;m<=3;m++){
for (n=0;n<=3;n++){
for (j=0;j<=3;j++){
G2[m][n] += -(Pijkl(V, x, dx, j, m, n, j));}}}

for (j=0;j<=3;j++){
TG2 += G2[j][j];}

for (j=0;j<=3;j++){
G2[j][j] = G2[j][j] - TG2/2;}

for (m=0;m<=3;m++){
for (n=0;n<=3;n++){
for (j=0;j<=3;j++){
for (k=0;k<=3;k++){
G3[m][n] += Pijk(V, x, dx, k, j, n) * Pijk(V, x, dx, j, m, k);}}}}

for (j=0;j<=3;j++){
TG3 += G3[j][j];}

for (j=0;j<=3;j++){
G3[j][j] = G3[j][j] - TG3/2;}

for (m=0;m<=3;m++){
for (n=0;n<=3;n++){
for (j=0;j<=3;j++){
for (k=0;k<=3;k++){
G4[m][n] += -(Pijk(V, x, dx, k, j, k) * Pijk(V, x, dx, j, m, n));}}}}

for (j=0;j<=3;j++){
TG4 += G4[j][j];}

for (j=0;j<=3;j++){
G4[j][j] = G4[j][j] - TG4/2;}


fputs("G1jk =",stdprn);fputs("\n",stdprn);
fputs("G1jk =",stdout);fputs("\n",stdout);
for(j=0;j<=3;j++){
for(k=0;k<=3;k++){
fputs(gcvt(G1[j][k], 12, buf), stdout); fputs(" ",stdout);
fputs(buf, stdprn); fputs(" ",stdprn);} fputs("\n",stdprn);fputs("\n",stdout);}

fputs("G2jk =",stdout);fputs("\n",stdout);
fputs("G2jk =",stdprn);fputs("\n",stdprn);
for(j=0;j<=3;j++){
for(k=0;k<=3;k++){
fputs(gcvt(G2[j][k], 12, buf), stdout); fputs(" ",stdout);
fputs(buf, stdprn); fputs(" ",stdprn);} fputs("\n",stdprn);fputs("\n",stdout);}

fputs("G3jk =",stdout);fputs("\n",stdout);
fputs("G3jk =",stdprn);fputs("\n",stdprn);
for(j=0;j<=3;j++){
for(k=0;k<=3;k++){
fputs(gcvt(G3[j][k], 12, buf), stdout); fputs(" ",stdout);
fputs(buf, stdprn); fputs(" ",stdprn);} fputs("\n",stdprn);fputs("\n",stdout);}

fputs("G4jk =",stdout);fputs("\n",stdout);
fputs("G4jk =",stdprn);fputs("\n",stdprn);
for(j=0;j<=3;j++){
for(k=0;k<=3;k++){
fputs(gcvt(G4[j][k], 12, buf), stdout); fputs(" ",stdout);
fputs(buf, stdprn); fputs(" ",stdprn);} fputs("\n",stdprn);fputs("\n",stdout);}
#else FOUR_COMPONENTS
fputs("FOUR COMPONENTS non defined", stdout);
#endif FOUR_COMPONENTS
goto a00;
a15:
;
}



Saturday the 19th. Joomla 2.5 Templates. Site built by Levi-Lidar: levilidar@gmail.com
Copyright 2012

©