/*R modified 28.5.91 for "Experimental ..." special calculations*/
#include <stdio.h>
#include <math.h>
#define DEF1
#define mrs (1 - rs/x[1])
#define mrs2 sqrt(mrs)
#define Z0 (double x[]){return 0;}
#define Z1 (double x[]){return 1;}
#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]) /*Lorentz parameters Boosts*/
#define B2 (Lp2/x[1])
#define B3 (Lp3/x[1])
#define B (Lp4/x[1]) /*Rotations*/
#define R2 (Lp5/x[1])
#define R3 (Lp6/x[1])
#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*/*/
#define AnaliticmixRM00 ((1/mrs2 - 1) * rs/(x[1]*x[1]*x[1]))
#define AnaliticmixRM11 0
#define AnaliticmixRM22 (0.5*(1/mrs2 - 1) * rs/(x[1]*x[1]*x[1])-(1-mrs2)*(1-mrs2
)/(x[1]*x[1]))
#define AnaliticmixRM33 AnaliticmixRM22
/*Added 28.5.91*/
#define AnaliticontRM00 mrs * ((1/mrs2 - 1) * rs/(x[1]*x[1]*x[1]))
#define AnaliticontEg0 (- rs/Gn/2 * (log(2*mrs2/(1+mrs2)) + 1/mrs2 - 1))
/*End Added*/
#define En (-0.5*Gn*M*M/r) /*Newton Gravitational Energy*/
#define Er (0.5*M*log(2*mrs/(0.5+0.5*mrs+sqrt(mrs)))) /*Relativistic integral gr
avitational energy kTg = Rg/2*//*Used symbols t,r,a,b,rs,mrs,mrs2,Z0,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); r
eturn z;}
#endif

d Lp1 = 0; /* Lorenz parameter - Boosts*/
d Lp2 = 0;
d Lp3 = 0;
d Lp4 = 0;/*Rotations*/
d Lp5 = 0;
d Lp6 = 0;

int ETA[4] = {1,-1,-1,-1};
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 original 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) * r
s/(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,A
naliticGG22,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}
;
/*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,F
0,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,F
0,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,L
B300};

/*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);}
/*K7.CAN was taken out from here*/

/*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(contd
erV,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;i
nt 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;}

/*Contravariant NGGlm*/
/* NGGlm = Val;i ETAaa Vai;m - Val;i ETAaa Vam;i */
double NGGlm(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,l,k) * ETA[j] * contderV(V,x,dx,j,k,m) * gg[m](x) -
contderV(V,x,dx,j,l,k) *ETA[j] * contderV(V,x,dx,j,m,k) * gg[k](x);}}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,d
x,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)*/

/*Added 28.5.91*/
double AncontRM00(double rs, double r){ double z; z=AnaliticontRM00; return z;}
double AncontEg0(double rs, double r){ double z; z=AnaliticontEg0; return z;}

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 NGG[4][4];
d L1 = 7e10;
d L2 = 1e12;
d idx = 1e-15;
d zz[4];
/*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*/
fputs("This version differs from Q in that GM was replaced by RM\n",stdout);
fputs(" Vai;jk was replaced by AnaliticmixRMjj/2\n",stdout);
fputs(" and 6 was replaced by NGGlm\n",stdout);
fputs("13 Integral is calculated with TG = 1/2k RG\n",stdout);
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);

/*K1.CAN was taken out from here*/
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("rs = ",stdout); fputs(gcvt(rs,12,buf),stdout);
fputs(" r = ",stdout); fputs(gcvt(r,12,buf),stdout);
/*K2.CAN was taken out from here*/

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);
/*K3.CAN was taken out from here*/

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:
/*K4.CAN was taken out from here*/

goto nl;/*Lorentz transformation disabled -Added 28.5.91*/
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 AnaliticmixRM00/2 6 NGG\n",stdout)
;
fputs("7 Rjk 8 RMjk 9 (VaiVbi;j);j 10 GGjk 11 AnaliticGGjk\n",stdout);
fputs(" 12 Contravariant TG00 13 Integral 14 FOUR COMPONENTS\n15 EXIT\n",stdo
ut);
gets(buf);
/*K5.CAN was taken out from here*/

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;
goto a00;
/*K6.CAN was taken out from here*/
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;
#ifdef DEF2
/* Fuera de Servicio*/
z = mrs;
fputs("mrs = ",stdout);
fputs(gcvt(z, 12, buf), stdout);
fputs("\n", stdout);
z = mrs2;
fputs("mrs2 = ", stdout);
fputs(gcvt(z, 12, buf), stdout);
fputs("\n",stdout);
fputs("g = \n", stdout);
for (j=0;j<=3;j++){
z = g[j](x);

fputs(gcvt(z, 12, buf), stdout);
fputs(" ", stdout);
}
fputs("\n\n", stdout);
#endif DEF2
/*Fin de "Fuera de Servicio"*/
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:

fputs(" j\n",stdout);
fputs("AnaliticmixRM /2 =\n",stdout);
fputs(" j\n",stdout);
fputs("AnaliticmixRMjj = \n",stdprn);
zz[0] = AnaliticmixRM00/2;zz[1]=AnaliticmixRM11;zz[2]=AnaliticmixRM22/2;zz[3]=A
naliticmixRM33/2;
for(j=0;j<=3;j++){
fputs(gcvt(zz[j],12,buf), stdout);fputs(" ",stdout);
fputs(buf, stdprn);fputs(" ",stdprn);
}
fputs("\n",stdout);
fputs("\n",stdprn);
goto a00;
/* 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:

for (j=0;j<=3;j++){
for (k=0;k<=3;k++){
NGG[j][k] = NGGlm(V,x,dx,j,k);
}}

fputs(" jk \n",stdout);
fputs("NGG =\n",stdout);
fputs("NGG =\n",stdprn);
fjk fputs(gcvt(NGG[j][k],12,buf),stdout); fputs(" ",stdout);
fputs(buf,stdprn);fputs(" ",stdprn);}
fputs("\n",stdout);
fputs("\n",stdprn);}
/* 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("RM =\n",stdout);
fputs(" k\n",stdout);
fputs("RM =\n",stdprn);
fjk fputs(gcvt(GM[j][k],12,buf),stdout); fputs(" ",stdout);
fputs(buf,stdprn);fputs(" ",stdprn);}
fputs("\n",stdout);
fputs("\n",stdprn);}

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
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("Analitical Calculation of Er when kTg = Rg/2) \n",stdout);
fputs("Er = sqrt(-g) V00 TG00 = ",stdout);
fputs("Er = sqrt(-g) V00 TG00 = ",stdprn);

fputs(gcvt(Er,12,buf), stdout);fputs(" ",stdout);
fputs(buf, stdprn);fputs(" ",stdprn);
fputs("\n",stdout);
fputs("\n",stdprn);
fputs("Er/En = ",stdout);fputs("E/En = ",stdprn);
fputs(gcvt(Er/En,12,buf), stdout);
fputs(gcvt(Er/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:

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

©