/*Analitic Calculations and Graphics of "Experimental ..."*/
/* Invariant and Covariant Integrals*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
/*#include <alloc.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)) + 2/mrs2 - 1/mrs -1)/*dEg/dM*/
/*#define V sqrt(1-(1+mrs*DEg)*(1+mrs*DEg))*/
#define V sqrt(1-(1+mrs*DEg)*(1+mrs*DEg)/mrs)

 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];
/*Initiate Graph*/
	DDc=0;
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();
/*Initiate line parameters*/
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;
/*Plotting*/
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();


/*Axis*/
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:

/*END GRAPH()*/



a1:
	fputs("FOUNDATI analitical calculations\n", stdout);


/*Ask for parameter values*/
/*	fputs("rs = ",stdout);
	 fputs(gcvt(rs, 12, buf),stdout);fputs("  rs = ",stdout);
	gets(buf);
	if (buf[0] != '')
	rs = atof(buf);

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

/*Print parameter values*/
/*	fputs("rs = ",stdout); fputs(gcvt(rs,12,buf),stdout);
	fputs("rs = ",stdprn); fputs(buf,stdprn);
	fputs("  R/rs = ",stdout); fputs(gcvt(R,12,buf),stdout);
        fputs("  R/rs = ",stdprn); fputs(buf,stdprn);
	fputs("\n",stdout);
	fputs("\n",stdprn);*/

a00:
fputs("\n1 Paramters  2 Velocity  3 v(dEg/dM) 4 Contravariant Eg/M 5 Invariant E
mg/M\n",stdout);
fputs("6 Covariant Mg/M  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);
		goto a00;
a3:
fputs("EMPTY\n", stdout);
		goto a00;

a4:
/*Conserved Gravitational Energy Integral*/
fputs("Eg/M = 2(ln(2 mrs2/(1+mrs2)) + 1/(1+mrs2) - 1/2)\n",stdout);
	fputs("Eg/M = ", stdout);fputs("Eg/M = ", stdprn);
	fputs(gcvt(EgM,12,buf),stdout);fputs(buf,stdprn);
	fputs("    En/M = ", stdout);fputs("    En/M = ", stdprn);
	fputs(gcvt(EnM,12,buf),stdout);fputs(buf,stdprn);
	fputs("\n", stdout);fputs("\n", stdprn);
	goto a1;
a5:
/* Invariant Gravitational Energy Integral*/
fputs("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("    En/M = ", stdout);fputs("    En/M = ", stdprn);
	fputs(gcvt(EnM,12,buf),stdout);fputs(buf,stdprn);
	fputs("\n", stdout);fputs("\n", stdprn);
	goto a1;
a6:
/*Covariant  Gravitational Energy Integral*/
fputs("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("    En/M = ", stdout);fputs("    En/M = ", stdprn);
	fputs(gcvt(EnM,12,buf),stdout);fputs(buf,stdprn);
	fputs("\n", stdout);fputs("\n", stdprn);
	goto a1;
E: