/*Analitic Calculations and Graphics of Experimental Facts and ..."*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "prngrph.h"
#define mrs (1.0 - 1.0/r)
#define mrs2 sqrt(mrs)
#define d double
#define rz return z;}
#define PI 3.1415927
#define K (8*PI*7.425*1e-29) /*GN/c2  Einstein Gravitational Constant in cm/gr*/
#define Gn 7.425e-29
#define M (1/(2*Gn))
#define R r /*Source Radious */
#define EnM (-0.25/R) /*Newton Gravitational Energy per unite mass*/
#define EgM (2*(log(2*mrs2/(1+mrs2)) + 1/(1+mrs2) - 0.5)) /* Conserved Gravitati
ona Energy per unite mass Integral*/
#define EmgM (-(1-mrs2)/(1+mrs2))/*Invariant Gravitational Energy per unite mass
 Integral*/
#define MgM  (2*(1/(1+mrs2)+log(1+mrs2)-0.5-log(2)))/*Covariant Gravitational En
ergy per unite mass Integral*/
#define Ve sqrt(1.0/r)
#define DEg (log(4*mrs/(1+mrs2)/(1+mrs2)) + (1-mrs)/(mrs2+mrs) + 1/mrs2 - 1/mrs)
/*d(contEg0)/dM)*/
#define V sqrt(1 - (1+mrs*DEg)*(1+mrs*DEg))
#define P (mrs2/(1+mrs*DEg))/*Scalar Dust Energy*/
/*DEg2 and V2 obtained from Tg=1/2 Rg*/
#define DEg2 (1 - log(2*mrs2/(1+mrs2)) - 1/mrs2 + (1-mrs2)/2/mrs - (1+mrs2)*(1-m
rs)/4/mrs/mrs)
#define V2 sqrt(1 - (1+mrs*DEg2)*(1+mrs*DEg2))/*From Tg=1/2 Rg*/


 void main (void)
	{
	int i;
	int j;
	double x;
	double y;
	double y1;
	double r=1;
	double z;
	double u;
	double v;
	double r0=0.0; double EmgM0=-2.0;

	double BeginX = 1.01;
	double EndX = 5.2;
	int full1 = 1;
	int doble1 = 0;
	int j1 = 1000;
	int full2 = 2;
	int doble2 = 0;
	int j2 = 50;
	int full3 = 0;
	int doble3 = 2;
	int j3 = 2000;
	int full4 = 1;
	int doble4 = 0;
	int j4 = 50;
	int full5 = 0;
	int doble5 = 0;
	int j5 = 50;

	static char buf[30];
np0:
	fputs("\n1 Do Graph 2 Get Graph 3 Analitical Calculations E  EXIT\n",stdout);
	gets(buf);
	if(buf[0]=='1') goto np1;
	if(buf[0]=='2') goto GD;
	if(buf[0]=='3') goto a1;
	if(buf[0]=='E') goto E;
	fputs("The only accepted characters are 1, 2, 3 or E",stdout); goto np0;

np1:
	run_initgraph();
np2:
	run_line_parameters();
np3:

fputs("\nGraph Menu\n",stdout);
fputs("\n1 Conserved Eg/M and En/M  2 Velocities \n",stdout);
fputs("3 Graph Parameters  4 Line Parameters  C Analitic Calculation E Exit\n",s
tdout); gets(buf);
	if(buf[0]=='1') goto l1;
	if(buf[0]=='2') goto l2;
	if(buf[0]=='3') {closegraph(); goto np1;}
	if(buf[0]=='4')  goto np2;
	if(buf[0]=='C') {closegraph(); goto a1;}
	if(buf[0]=='E') {closegraph(); goto E;}
fputs("\nDigits 1 to 4, and C and E are the only accepted entries\n",stdout);
	goto np3;
l1:
/*Conserved Relativistic Gravitational Energy for Unit Mass */
fputs("\nGraph of the Conserved Gravitational Energy Eg/M and Newton Gravitation
al Energy\n",stdout);
fputs("En/M per Unit Mass against R/2MG\n",stdout);
fputs("\nLine Parameters for Eg/M (Contravariant)\n",stdout);
fputs("Default/Previous Values\n",stdout);

fputs("Begin with R/2MG = ",stdout);fputs(gcvt(x=BeginX,10,buf),stdout);
fputs("  End with R/2MG = ",stdout); fputs(gcvt(y1=EndX,10,buf),stdout);
fputs("\nfull = ",stdout); fputs(itoa(full=full1,buf,10),stdout);
fputs("   doble = ",stdout); fputs(itoa(doble=doble1,buf,10),stdout);
fputs("   Number of Points = ",stdout); fputs(itoa(j=j1,buf,10),stdout);
fputs("\nOkey?(y/n)\n",stdout);gets(buf);
	if(buf[0]=='y') goto l11;
fputs("<ENTER> to confirm previous value\n",stdout);
again1:
fputs("Begin with R/2MG = ",stdout);gets(buf);if(buf[0]!='')BeginX=x=atof(buf);
if(x<minX || x>maxX) {fputs("R/2MG has to be between minX and maxX\n",stdout);
			goto again1;}
again2:
fputs("End with R/2MG = ",stdout);gets(buf);if(buf[0]!='')EndX=y1=atof(buf);
if(y1<minX || y1>maxX) {fputs("R/2MG has to be between minX and maxX\n",stdout);
			goto again2;}
fputs("full = (0,1 or 2) ",stdout);gets(buf); if(buf[0]!='')full=full1=atoi(buf)
;
fputs("doble = (0,1 or 2) ",stdout);gets(buf); if(buf[0]!='');doble=doble1=atoi(
buf);
fputs("Number of points = ",stdout); gets(buf);if(buf[0]!='')j=j1=atoi(buf);
l11:
fputs("\nShould I print first line parameters (Eg/M)?(y/n)",stdout);gets(buf);
	if(buf[0]=='n') goto l12;
	fputs("\nFirst Line Parameters (Eg/M)\n",stdprn);
	fputs("Begin with R/2MG = ",stdprn);fputs(gcvt(BeginX,10,buf),stdprn);
	fputs("End with R/2MG = ",stdprn);fputs(gcvt(EndX,10,buf),stdprn);
	fputs("Number of points = ",stdprn); fputs(itoa(j,buf,10),stdprn);
	fputs("   full = ",stdprn); fputs(itoa(full,buf,10),stdprn);
	fputs("   doble = ",stdprn); fputs(itoa(doble,buf,10),stdprn);fputs("\n",stdprn
);
l12:
	for(i=0;i<j;i++){
	r = x + (y1-1)/(j-1)*i;
	recpoint(r, EgM);}
	newline();

fputs("READING Eg/M line characteristic",stdout);
fputs("\ny = -1,75  (<ENTER> to comfirm)  y = ",stdout); y = -1.75; gets(buf);
	if(buf[0]!='') y = atof(buf);
fputs("Initial x = 3.0 (<ENTER> to comfirm) Initial x = ",stdout); x = 3.0; gets
(buf);
	if(buf[0]!='') x = atof(buf);
fputs("Final x = 4.0 (<ENTER> to comfirm) Final x = ",stdout); r = 4.0; gets(buf
);
	if(buf[0]!='') r = atof(buf);
fputs("Number of Points = 10 (<ENTER> to comfirm) No. of Points = ",stdout); j =
 10; gets(buf);
	if(buf[0]!='') j = atoi(buf);
fputs("\nShould I print READING first line parameters?(y/n)",stdout);gets(buf);
	if(buf[0]=='n') goto l19;
	fputs("\nREADING First Line Parameters\n",stdprn);
	fputs("Number of points = ",stdprn); fputs(itoa(j,buf,10),stdprn);
	fputs("   y = ",stdprn); fputs(gcvt(y,10,buf),stdprn);
	fputs("   Initial x = ",stdprn); fputs(gcvt(x,10,buf),stdprn);
	fputs("   Final x = ",stdprn); fputs(gcvt(r,10,buf),stdprn);fputs("\n",stdprn);
l19:
	for(i=0;i<j;i++){
	z = x + (r-x)/(j-1)*i;
	recpoint(z, y);}
	newline();
/*ARROW*/
	fputs("\nArrow x = 1.025 (<ENTER> to confirm)  x = ",stdout); x = 1.025;
	gets(buf);if(buf[0]!='') x = atof(buf);
	fputs("\nArrow y = -1.97 (<ENTER> to confirm)  y = ",stdout); y = -1.97;
	gets(buf);if(buf[0]!='') y = atof(buf);
	full = 1;
	recpoint(x, y);
	recpoint(x, minY);
	recpoint(1, y);
	recpoint(x+x-1, y);
	recpoint(x, minY);
	newline();
/*Second line*/
fputs("\nLine Parameters for En/M (Invariant)\n",stdout);
fputs("Default/Previous Values\n",stdout);
fputs("full = ",stdout); fputs(itoa(full=full2,buf,10),stdout);
fputs("   doble = ",stdout); fputs(itoa(doble=doble2,buf,10),stdout);
fputs("   Number of Points = ",stdout); fputs(itoa(j=j2,buf,10),stdout);
fputs("\nOkey?(y/n)\n",stdout);gets(buf);
	if(buf[0]=='y') goto l21;
fputs("<ENTER> to confirm previous value\n",stdout);
fputs("full = (0,1 or 2) ",stdout);gets(buf); if(buf[0]!='')full=full2=atoi(buf)
;
fputs("doble = (0,1 or 2) ",stdout);gets(buf); if(buf[0]!='');doble=doble2=atoi(
buf);
fputs("Number of points = ",stdout); gets(buf);if(buf[0]!='')j=j2=atoi(buf);
l21:
fputs("\nShould I print second line parameters (En/M)?(y/n)",stdout);gets(buf);
	if(buf[0]=='n') goto l22;
	fputs("\nSecond Line Parameters\n",stdprn);
	fputs("Number of points = ",stdprn); fputs(itoa(j,buf,10),stdprn);
	fputs("   full = ",stdprn); fputs(itoa(full,buf,10),stdprn);
	fputs("   doble = ",stdprn); fputs(itoa(doble,buf,10),stdprn);fputs("\n",stdprn
);
l22:
	for(i=0;i<j;i++){
	r = 1 + (y1-1)/(j-1)*i;
	recpoint(r, EnM);}
	newline();

fputs("READING En/M line characteristic",stdout);
fputs("\ny = -1.03  (<ENTER> to comfirm)  y = ",stdout); y = -1.00; gets(buf);
	if(buf[0]!='') y = atof(buf);
fputs("Initial x = 3.0 (<ENTER> to comfirm) Initial x = ",stdout); x = 3.0; gets
(buf);
	if(buf[0]!='') x = atof(buf);
fputs("Final x = 4.0 (<ENTER> to comfirm) Final x = ",stdout); r = 4.0; gets(buf
);
	if(buf[0]!='') r = atof(buf);
fputs("Number of Points (<ENTER> to comfirm) = ",stdout);
fputs(itoa(j=j/4,buf,10),stdout);
fputs("No. of Points = ",stdout); gets(buf);if(buf[0]!='') j = atoi(buf);
fputs("\nShould I print Reading second line parameters?(y/n)",stdout);gets(buf);
	if(buf[0]=='n') goto l201;
	fputs("\nREADING Second Line Parameters\n",stdprn);
	fputs("Number of points = ",stdprn); fputs(itoa(j,buf,10),stdprn);
	fputs("   y = ",stdprn); fputs(gcvt(y,10,buf),stdprn);
	fputs("   Initial x = ",stdprn); fputs(gcvt(x,10,buf),stdprn);
	fputs("   Final x = ",stdprn); fputs(gcvt(r,10,buf),stdprn);fputs("\n",stdprn);
l201:
	for(i=0;i<j;i++){
	z = x + (r-x)/(j-1)*i;
	recpoint(z, y);}
	newline();
/*Third line*/
fputs("\nLine Parameters for Emg/M (Invariant)\n",stdout);
fputs("Default/Previous Values\n",stdout);
fputs("full = ",stdout); fputs(itoa(full=full3,buf,10),stdout);
fputs("   doble = ",stdout); fputs(itoa(doble=doble3,buf,10),stdout);
fputs("   Number of Points = ",stdout); fputs(itoa(j=j3,buf,10),stdout);
fputs("\nOkey?(y/n)\n",stdout);gets(buf);
	if(buf[0]=='y') goto l31;
fputs("<ENTER> to confirm previous value\n",stdout);
fputs("full = (0,1 or 2) ",stdout);gets(buf); if(buf[0]!='')full=full3=atoi(buf)
;
fputs("doble = (0,1 or 2) ",stdout);gets(buf); if(buf[0]!='');doble=doble3=atoi(
buf);
fputs("Number of points = ",stdout); gets(buf);if(buf[0]!='')j=j3=atoi(buf);
l31:
fputs("\nShould I print third line parameters (Emg)?(y/n)",stdout);gets(buf);
	if(buf[0]=='n') goto l32;
	fputs("\nThird Line Parameters (Invariant Emg/M)\n",stdprn);
	fputs("Number of points = ",stdprn); fputs(itoa(j,buf,10),stdprn);
	fputs("   full = ",stdprn); fputs(itoa(full,buf,10),stdprn);
	fputs("   doble = ",stdprn); fputs(itoa(doble,buf,10),stdprn);fputs("\n",stdprn
);
l32:
	r0=0.0;EmgM0=-2.0;
	for(i=0;i<j;i++){
	r = 1 + (y1-1)/(j-1)*i;
	if(r-r0>0.15||EmgM-EmgM0>0.07){r0=r;EmgM0=EmgM;
	recpoint(r, EmgM);}}
	newline();

fputs("Reading Emg/M line characteristic",stdout);
fputs("\ny = -1.51  (<ENTER> to comfirm)  y = ",stdout); y = -1.5; gets(buf);
	if(buf[0]!='') y = atof(buf);
fputs("Initial x = 3.0 (<ENTER> to comfirm) Initial x = ",stdout); x = 3.0; gets
(buf);
	if(buf[0]!='') x = atof(buf);
fputs("Final x = 4.0 (<ENTER> to comfirm) Final x = ",stdout); r = 4.0; gets(buf
);
	if(buf[0]!='') r = atof(buf);
fputs("Number of Points (<ENTER> to comfirm) = ",stdout);
fputs(itoa(j=j/4,buf,10),stdout);
fputs("No. of Points = ",stdout); gets(buf);if(buf[0]!='') j = atoi(buf);
fputs("\nShould I print Reading third line parameters?(y/n)",stdout);gets(buf);
	if(buf[0]=='n') goto l301;
	fputs("\nREADING Third Line Parameters Emg/M\n",stdprn);
	fputs("Number of points = ",stdprn); fputs(itoa(j,buf,10),stdprn);
	fputs("   y = ",stdprn); fputs(gcvt(y,10,buf),stdprn);
	fputs("   Initial x = ",stdprn); fputs(gcvt(x,10,buf),stdprn);
	fputs("   Final x = ",stdprn); fputs(gcvt(r,10,buf),stdprn);fputs("\n",stdprn);
l301:
	 r0=0.0;
	for(i=0;i<j;i++){
        z = x + (r-x)/(j-1)*i;
	if(z-r0>0.15){r0=z;
	recpoint(z, y);}}
	newline();

/*Fourth line*/
fputs("\nLine Parameters for Mg/M (Invariant)\n",stdout);
fputs("Default/Previous Values\n",stdout);
fputs("full = ",stdout); fputs(itoa(full=full4,buf,10),stdout);
fputs("   doble = ",stdout); fputs(itoa(doble=doble4,buf,10),stdout);
fputs("   Number of Points = ",stdout); fputs(itoa(j=j4,buf,10),stdout);
fputs("\nOkey?(y/n)\n",stdout);gets(buf);
	if(buf[0]=='y') goto l41;
fputs("<ENTER> to confirm previous value\n",stdout);
fputs("full = (0,1 or 2) ",stdout);gets(buf); if(buf[0]!='')full=full4=atoi(buf)
;
fputs("doble = (0,1 or 2) ",stdout);gets(buf); if(buf[0]!='');doble=doble4=atoi(
buf);
fputs("Number of points = ",stdout); gets(buf);if(buf[0]!='')j=j4=atoi(buf);
l41:
fputs("\nShould I print fourth line parameters (Mg/M)?(y/n)",stdout);gets(buf);
	if(buf[0]=='n') goto l42;
	fputs("\nFourth Line Parameters (Mg/M)\n",stdprn);
	fputs("Number of points = ",stdprn); fputs(itoa(j,buf,10),stdprn);
	fputs("   full = ",stdprn); fputs(itoa(full,buf,10),stdprn);
	fputs("   doble = ",stdprn); fputs(itoa(doble,buf,10),stdprn);fputs("\n",stdprn
);
l42:
	for(i=0;i<j;i++){
	r = 1 + (y1-1)/(j-1)*i;
	if(i%3==2){newline();continue;}
	recpoint(r, MgM);}
	newline();

fputs("Reading Mg/M line characteristic",stdout);
fputs("\ny = -1.27  (<ENTER> to comfirm)  y = ",stdout); y = -1.25; gets(buf);
	if(buf[0]!='') y = atof(buf);
fputs("Initial x = 3.0 (<ENTER> to comfirm) Initial x = ",stdout); x = 3.0; gets
(buf);
	if(buf[0]!='') x = atof(buf);
fputs("Final x = 4.0 (<ENTER> to comfirm) Final x = ",stdout); r = 4.0; gets(buf
);
	if(buf[0]!='') r = atof(buf);
fputs("Number of Points (<ENTER> to comfirm) = ",stdout);
fputs(itoa(j=j/4,buf,10),stdout);
fputs("No. of Points = ",stdout); gets(buf);if(buf[0]!='') j = atoi(buf);
fputs("\nShould I print Reading fourth (MgM) line parameters?(y/n)",stdout);gets
(buf);
	if(buf[0]=='n') goto l401;
	fputs("\nREADING Fourth Line Parameters Mg/M\n",stdprn);
	fputs("Number of points = ",stdprn); fputs(itoa(j,buf,10),stdprn);
	fputs("   y = ",stdprn); fputs(gcvt(y,10,buf),stdprn);
	fputs("   Initial x = ",stdprn); fputs(gcvt(x,10,buf),stdprn);
	fputs("   Final x = ",stdprn); fputs(gcvt(r,10,buf),stdprn);fputs("\n",stdprn);
l401:
	for(i=0;i<j;i++){
	z = x + (r-x)/(j-1)*i;
	if(i%3==2){newline();continue;}
	recpoint(z, y);}
	newline();

/*Fifth line - x = 1*/
fputs("Line Parameters for the vertical x = 1  - <ENTER> to confirm",stdout);
fputs("\nfull = 0 ",stdout); full = 0; gets(buf);if(buf[0]!='')full=atoi(buf);
fputs("doble = 0  doble = ",stdout); doble = 0;
gets(buf); if(buf[0]!='') doble=atoi(buf);
fputs("Number of points = 50 (<ENTER> to confirm)  N. of Points = ",stdout);
	j=50; gets(buf);if(buf[0]!='') j = atoi(buf);
fputs("\nShould I print fifth line parameters?(y/n)",stdout);gets(buf);
	if(buf[0]=='n') goto l17;
	fputs("\nFifth Line Parameters (Vertical x = 1\n",stdprn);
	fputs("Number of points = ",stdprn); fputs(itoa(j,buf,10),stdprn);
	fputs("   full = ",stdprn); fputs(itoa(full,buf,10),stdprn);
	fputs("   doble = ",stdprn); fputs(itoa(doble,buf,10),stdprn);fputs("\n",stdprn
);
l17:
	for(i=0;i<=j;i++){
	r = minY/j*i;
	recpoint(1, r);}
	newline();
Axis:
fputs("\nAxis Parameters\n",stdout);
fputs("\nAxis?(y/n)",stdout);gets(buf);
	if(buf[0]=='n') goto l13;
	fputs("Y axis coordinate = ",stdout);fputs(gcvt(x=0.0,10,buf),stdout);
	fputs("  Y axis coordinate = ",stdout);gets(buf);
	if(buf[0]!='')x = atof(buf);
	fputs("X axis abscissa = ",stdout);fputs(gcvt(y=0.0,10,buf),stdout);
	fputs("  X axis abscissa = ",stdout);gets(buf);
	if(buf[0]!='')y = atof(buf);
	fputs("Xmark = ",stdout); fputs(gcvt(u=1.0,10,buf),stdout);
	fputs("  Xmark = ",stdout);gets(buf);
	if(buf[0]!='')u = atof(buf);
	fputs("Ymark = ",stdout);fputs(gcvt(v=0.25,10,buf),stdout);
	fputs("  Ymark = ",stdout);gets(buf);
	if(buf[0]!='')v = atof(buf);
/*Print Axes Parameters*/
	fputs("Yaxis = ",stdout); fputs(gcvt(x,10,buf),stdout);
	fputs("   Xaxis = ",stdout); fputs(gcvt(y,10,buf),stdout);
	fputs("   Xmark = ",stdout); fputs(gcvt(u,10,buf),stdout);
	fputs("   Ymark = ",stdout); fputs(gcvt(v,10,buf),stdout);
	fputs("\nShould I print the axes parameters?(y/n)",stdout);gets(buf);
	if(buf[0]=='n') goto l15;
	fputs("\nYaxis = ",stdprn); fputs(gcvt(x,10,buf),stdprn);
	fputs("   Xaxis = ",stdprn); fputs(gcvt(y,10,buf),stdprn);
	fputs("   Xmark = ",stdprn); fputs(gcvt(u,10,buf),stdprn);
	fputs("   Ymark = ",stdprn); fputs(gcvt(v,10,buf),stdprn);fputs("\n",stdprn);
l15:
	axis(x, y, u, v, 0);
	newline();

l13:
	redata();/*Reduces the data*/
fputs("\nEnter Filename to Save or nn\n",stdout); gets(buf);
	if(buf[0]=='n'&&buf[1]=='n') goto l16;
	savefile =(FILE *) buf;
	savegraph(savefile);
l16:
	printgraph();
	closegraph();
	goto np0;
GD:
fputs("\nEnter Filename to Print or nn\n",stdout); gets(buf);
	if(buf[0]=='n'&&buf[1]=='n') goto np0;
	savefile = (FILE *) buf;
	getgraph(savefile);
        printgraph();
	closegraph();
	goto np0;
l2:
/*Graph of compared velocities*/
fputs("\nGraph of Compared Velocities\n",stdout);
fputs("V/c against R/2MG of a Falling Body\n",stdout);
fputs("\nLine Parameters for Ve/C (Einstein Velocity)\n",stdout);
fputs("Default/Previous Values\n",stdout);

fputs("maxX = ",stdout);fputs(gcvt(maxX = 5.2, 10, buf),stdout);
fputs("  minX = ",stdout);fputs(gcvt(minX = -0.05, 10, buf),stdout);
fputs("  maxY = ",stdout);fputs(gcvt(maxY = 1.05, 10, buf),stdout);
fputs("  minY = ",stdout);fputs(gcvt(minY = -0.05, 10, buf),stdout);
fputs("\nBegin with R/2MG = ",stdout);fputs(gcvt(x=BeginX = 1, 10, buf),stdout);
fputs("  End with R/2MG = ",stdout); fputs(gcvt(y1=EndX=5.2 ,10,buf),stdout);
fputs("\nfull = ",stdout); fputs(itoa(full=full1,buf,10),stdout);
fputs("   doble = ",stdout); fputs(itoa(doble=doble1,buf,10),stdout);
fputs("   Number of Points = ",stdout); fputs(itoa(j=50,buf,10),stdout);
fputs("\nOkey?(y/n)\n",stdout);gets(buf);
	if(buf[0]=='y') goto l11x;
fputs("<ENTER> to confirm previous value\n",stdout);
fputs("maxX = ",stdout);gets(buf);if(buf[0]!='')maxX=atof(buf);
fputs("minX = ",stdout);gets(buf);if(buf[0]!='')minX=atof(buf);
fputs("maxY = ",stdout);gets(buf);if(buf[0]!='')maxY=atof(buf);
fputs("minY = ",stdout);gets(buf);if(buf[0]!='')minY=atof(buf);

again11:
fputs("Begin with R/2MG = ",stdout);gets(buf);if(buf[0]!='')BeginX=x=atof(buf);
if(x<minX || x>maxX) {fputs("R/2MG has to be between minX and maxX\n",stdout);
			goto again11;}
again22:
fputs("End with R/2MG = ",stdout);gets(buf);if(buf[0]!='')EndX=y1=atof(buf);
if(y1<minX || y1>maxX) {fputs("R/2MG has to be between minX and maxX\n",stdout);
			goto again22;}
fputs("full = (0,1 or 2) ",stdout);gets(buf); if(buf[0]!='')full=full1=atoi(buf)
;
fputs("doble = (0,1 or 2) ",stdout);gets(buf); if(buf[0]!='');doble=doble1=atoi(
buf);
fputs("Number of points = ",stdout); gets(buf);if(buf[0]!='')j=atoi(buf);
l11x:
fputs("\nShould I print Geodesic Path Velocity line parameters (Ve/c)?(y/n)",std
out);gets(buf);
	if(buf[0]=='n') goto l12x;
	fputs("\nGeodesic Path Velocity Line Parameters (Ve/c)\n",stdprn);
	fputs("maxX = ",stdprn);fputs(gcvt(maxX, 10, buf),stdprn);
fputs("  minX = ",stdprn);fputs(gcvt(minX, 10, buf),stdprn);
fputs("  maxY = ",stdprn);fputs(gcvt(maxY, 10, buf),stdprn);
fputs("  minY = ",stdprn);fputs(gcvt(minY, 10, buf),stdprn);
	fputs("Begin with R/2MG = ",stdprn);fputs(gcvt(BeginX,10,buf),stdprn);
	fputs("End with R/2MG = ",stdprn);fputs(gcvt(EndX,10,buf),stdprn);
	fputs("Number of points = ",stdprn); fputs(itoa(j,buf,10),stdprn);
	fputs("   full = ",stdprn); fputs(itoa(full,buf,10),stdprn);
	fputs("   doble = ",stdprn); fputs(itoa(doble,buf,10),stdprn);fputs("\n",stdprn
);
l12x:
	for(i=0;i<j;i++){
        r = x + (y1-1)/(j-1)*i;
	if(i%3==2){newline();continue;}
	recpoint(r, Ve);}
	newline();

fputs("READING G.Path Velocity line characteristic",stdout);
fputs("\ny = 0.60  (<ENTER> to comfirm)  y = ",stdout); y = 0.60; gets(buf);
	if(buf[0]!='') y = atof(buf);
fputs("Initial x = 3.0 (<ENTER> to comfirm) Initial x = ",stdout); x = 3.0; gets
(buf);
	if(buf[0]!='') x = atof(buf);
fputs("Final x = 4.0 (<ENTER> to comfirm) Final x = ",stdout); r = 4.0; gets(buf
);
	if(buf[0]!='') r = atof(buf);
fputs("Number of Points = 10 (<ENTER> to comfirm) No. of Points = ",stdout); j =
 10; gets(buf);
	if(buf[0]!='') j = atoi(buf);
fputs("\nShould I print READING G.P.Velocity line parameters?(y/n)",stdout);gets
(buf);
	if(buf[0]=='n') goto l19x;
	fputs("\nREADING G.P.Velocity Line Parameters\n",stdprn);
	fputs("Number of points = ",stdprn); fputs(itoa(j,buf,10),stdprn);
	fputs("   y = ",stdprn); fputs(gcvt(y,10,buf),stdprn);
	fputs("   Initial x = ",stdprn); fputs(gcvt(x,10,buf),stdprn);
	fputs("   Final x = ",stdprn); fputs(gcvt(r,10,buf),stdprn);fputs("\n",stdprn);
l19x:
	for(i=0;i<j;i++){
	z = x + (r-x)/(j-1)*i;
	if(i%3==2){newline();continue;}
	recpoint(z, y);}
	newline();

/*Energy Conservation Velocity*/
fputs("\nLine Parameters for V/C (Energy Conservation Velocity)\n",stdout);
fputs("Default/Previous Values\n",stdout);
fputs("maxX = ",stdout);fputs(gcvt(maxX, 10, buf),stdout);
fputs("  minX = ",stdout);fputs(gcvt(minX, 10, buf),stdout);
fputs("  maxY = ",stdout);fputs(gcvt(maxY, 10, buf),stdout);
fputs("  minY = ",stdout);fputs(gcvt(minY, 10, buf),stdout);
fputs("\nBegin with R/2MG = ",stdout);fputs(gcvt(x=BeginX, 10, buf),stdout);
fputs("  End with R/2MG = ",stdout); fputs(gcvt(y1=EndX,10,buf),stdout);
fputs("\nfull = ",stdout); fputs(itoa(full=full1,buf,10),stdout);
fputs("   doble = ",stdout); fputs(itoa(doble=doble1,buf,10),stdout);
fputs("   Number of Points = ",stdout); fputs(itoa(j=50,buf,10),stdout);
fputs("\nOkey?(y/n)\n",stdout);gets(buf);
	if(buf[0]=='y') goto l11y;
fputs("<ENTER> to confirm previous value\n",stdout);
fputs("maxX = ",stdout);gets(buf);if(buf[0]!='')maxX=atof(buf);
fputs("minX = ",stdout);gets(buf);if(buf[0]!='')minX=atof(buf);
fputs("maxY = ",stdout);gets(buf);if(buf[0]!='')maxY=atof(buf);
fputs("minY = ",stdout);gets(buf);if(buf[0]!='')minY=atof(buf);

again1x:
fputs("Begin with R/2MG = ",stdout);gets(buf);if(buf[0]!='')BeginX=x=atof(buf);
if(x<minX || x>maxX) {fputs("R/2MG has to be between minX and maxX\n",stdout);
			goto again1x;}
again2x:
fputs("End with R/2MG = ",stdout);gets(buf);if(buf[0]!='')EndX=y1=atof(buf);
if(y1<minX || y1>maxX) {fputs("R/2MG has to be between minX and maxX\n",stdout);
			goto again2x;}
fputs("full = (0,1 or 2) ",stdout);gets(buf); if(buf[0]!='')full=full1=atoi(buf)
;
fputs("doble = (0,1 or 2) ",stdout);gets(buf); if(buf[0]!='');doble=doble1=atoi(
buf);
fputs("Number of points = ",stdout); gets(buf);if(buf[0]!='')j=atoi(buf);
l11y:
fputs("\nShould I print Energy Conservation Velocity line parameters (V/c)?(y/n)
",stdout);gets(buf);
	if(buf[0]=='n') goto l12y;
	fputs("\nEnergy Conservation velocity Line Parameters (Ve/c)\n",stdprn);
	fputs("maxX = ",stdprn);fputs(gcvt(maxX, 10, buf),stdprn);
fputs("  minX = ",stdprn);fputs(gcvt(minX, 10, buf),stdprn);
fputs("  maxY = ",stdprn);fputs(gcvt(maxY, 10, buf),stdprn);
fputs("  minY = ",stdprn);fputs(gcvt(minY, 10, buf),stdprn);
	fputs("Begin with R/2MG = ",stdprn);fputs(gcvt(BeginX,10,buf),stdprn);
	fputs("End with R/2MG = ",stdprn);fputs(gcvt(EndX,10,buf),stdprn);
	fputs("Number of points = ",stdprn); fputs(itoa(j,buf,10),stdprn);
	fputs("   full = ",stdprn); fputs(itoa(full,buf,10),stdprn);
	fputs("   doble = ",stdprn); fputs(itoa(doble,buf,10),stdprn);fputs("\n",stdprn
);
l12y:
	for(i=0;i<j;i++){
	r = x + (y1-x)/(j-1)*i;
	recpoint(r, V);}
	newline();

fputs("READING Energy Conservation Velocity line characteristic",stdout);
fputs("\ny = 0.50  (<ENTER> to comfirm)  y = ",stdout); y = 0.50; gets(buf);
	if(buf[0]!='') y = atof(buf);
fputs("Initial x = 3.0 (<ENTER> to comfirm) Initial x = ",stdout); x = 3.0; gets
(buf);
	if(buf[0]!='') x = atof(buf);
fputs("Final x = 4.0 (<ENTER> to comfirm) Final x = ",stdout); r = 4.0; gets(buf
);
	if(buf[0]!='') r = atof(buf);
fputs("Number of Points = 10 (<ENTER> to comfirm) No. of Points = ",stdout); j =
 10; gets(buf);
	if(buf[0]!='') j = atoi(buf);
fputs("\nShould I print READING E.C. Velocity line parameters?(y/n)",stdout);get
s(buf);
	if(buf[0]=='n') goto l19y;
	fputs("\nREADING E.C.Velocity Line Parameters\n",stdprn);
	fputs("Number of points = ",stdprn); fputs(itoa(j,buf,10),stdprn);
	fputs("   y = ",stdprn); fputs(gcvt(y,10,buf),stdprn);
	fputs("   Initial x = ",stdprn); fputs(gcvt(x,10,buf),stdprn);
	fputs("   Final x = ",stdprn); fputs(gcvt(r,10,buf),stdprn);fputs("\n",stdprn);
l19y:
	for(i=0;i<j;i++){
	z = x + (r-x)/(j-1)*i;
	recpoint(z, y);}
	newline();
	goto Axis;

a1:
	fputs("EMPTY\n", stdout);

a00:
fputs("\n1 EMPTY  2 Velocity  3 v(dEg/dM) 4 En/M Mg/M Emg/M Eg/M 5 P/dm\n",stdou
t);
fputs("6 EMPTY  E EXIT\n",stdout);
	gets(buf);

	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] == 'E') goto E;
	goto a00;

a2:
/*Velocities*/
	fputs("\nr/rs =  ",stdout);
	fputs(gcvt(r,12,buf),stdout);
	fputs("   r/rs =  ",stdout); gets(buf);if(buf[0]!='')r=atof(buf);
	fputs("\nr/rs =  ",stdprn);
	fputs(buf,stdprn);
	fputs("\nE.C. velocity V = ",stdout);fputs(gcvt(V,10,buf),stdout);
	fputs("  E.C. velocity V = ",stdprn);fputs(buf,stdprn);
	fputs(" G.P. velocity Ve = ",stdout);fputs(gcvt(Ve,10,buf),stdout);
	fputs("  G.P. velocity Ve = ",stdprn);fputs(buf,stdprn);fputs("\n",stdprn);
	fputs("\n",stdout);
	fputs("\nE.C. velocity from T=1/2K Rg  V2 = ",stdout);fputs(gcvt(V2,10,buf),std
out);
	fputs("  E.C. velocity from T=1/2K Rg V2 = ",stdprn);fputs(buf,stdprn);
	fputs("\n",stdout); fputs("\n",stdprn);
		goto a00;
a3:
fputs("EMPTY\n", stdout);
		goto a00;

a4:
/*Gravitational Energy Integrals*/
	fputs("\nr/rs =  ",stdout);
	fputs(gcvt(r,12,buf),stdout);
	fputs("   r/rs =  ",stdout); gets(buf);if(buf[0]!='')r=atof(buf);
	fputs("\nr/rs =  ",stdprn);
	fputs(buf,stdprn);
	fputs("    En/M = ", stdout);fputs("    En/M = ", stdprn);
	fputs(gcvt(EnM,12,buf),stdout);fputs(buf,stdprn);
fputs("\n(Covariant) Mg/M = 2 (1/1+mrs2 + ln(1+mrs2) - 0.5 - ln 2) \n",stdout);
	fputs("Mg/M = ", stdout);fputs("   Mg/M = ", stdprn);
	fputs(gcvt(MgM,12,buf),stdout);fputs(buf,stdprn);
        fputs("\n", stdout);fputs("\n", stdprn);
fputs("(Inv) Emg/M = - (1 - mrs2)/(1 + mrs2)\n",stdout);
	fputs("Emg/M = ", stdout);fputs("Emg/M = ", stdprn);
	fputs(gcvt(EmgM,12,buf),stdout);fputs(buf,stdprn);
fputs("\nConEg/M = 2(ln(2 mrs2/(1+mrs2)) + 1/(1+mrs2) - 1/2)\n",stdout);
	fputs("ConEg/M = ", stdout);fputs("   ConEg/M = ", stdprn);
	fputs(gcvt(EgM,12,buf),stdout);fputs(buf,stdprn);
	fputs("\n", stdout);fputs("\n", stdprn);

	goto a00;
a5:
	fputs("Scalar Energy of the Dust for Unite Rest Mass at Infinite P/dm", stdout)
;
	fputs("\n", stdout);
        fputs("\nr/rs =  ",stdout);
	fputs(gcvt(r,12,buf),stdout);
	fputs("   r/rs =  ",stdout); gets(buf);if(buf[0]!='')r=atof(buf);
	fputs("\nr/rs =  ",stdprn);
	fputs(buf,stdprn);
	fputs("    P/dm = ", stdout);fputs("    P/dm = ", stdprn);
	fputs(gcvt(P,12,buf),stdout);fputs(buf,stdprn);
        fputs("\n", stdout);fputs("\n", stdprn);
	goto a00;
a6:
	fputs("EMPTY", stdout);
	fputs("\n", stdout);
	goto a00;
E: