/* 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: ; }
כתיבת תגובה