/* Change Log
	5/10/86 	v2.2 Added single character aliases for satellite
			names.

	4/30/86 	v2.1 Print blank line if satellite dips below
			horizon and reappears during same orbit and day

	4/29/86 	v2.0 Changed GetSatelliteParams() to use AMsAT's
			"kepler.dat" file. Moved schedule to "mode.dat" file.

	4/22/86 	v1.3  Inserted N8FJB's suggestions for variable naming
			which maintain 8 character uniqueness.
			Also removed "include" file orbitr.h, which had two
			definitions of external functions defined in orbit.c
			    -K3MC 

	4/1/86		v1.2  Corrected a scanf conversion to %d for an int
			type.	 -K3MC
 
	3/19/86 	v1.1  Changed GetSatelliteParams to not pass NULL
			to sscanf.
 
									*/
 
#include "df2:include/lattice/stdio.h"
#include "df2:include/lattice/math.h"
 
extern double Kepler();
extern double GetDay();
 
 
/*   #define PI 3.14159265  */

#define PI2 (PI*2)
#define MinutesPerDay (24*60.0)
#define SecondsPerDay (60*MinutesPerDay)
#define HalfSecond (0.5/SecondsPerDay)
#define EarthRadius 6378.16		/* Kilometers		*/
#define C 2.997925e5			/* Kilometers/Second	*/
#define DegreesPerRadian (180/PI)
#define RadiansPerDegree (PI/180)
#define ABS(x) ((x) < 0 ? (-(x)) : (x))
#define SQR(x) ((x)*(x))
 
#define MaxModes 10
typedef struct {
		int MinPhase,MaxPhase;
		char ModeStr[20];
	       }  ModeRec;
  
static char Str1[30] ={ "N3EMO Orbit Simulator  v2.2"};
static char Str2[30] = {"N3DPU Version for Amiga "};
   

   /*	Keplerian Elements and misc. data for the satellite		 */

    double  EpochDay;			/* time of epoch		 */
    double  EpochMeanAnomaly;	/* Mean Anomaly at epoch	 */
    long    EpochOrbitNum; 		/* Integer orbit # of epoch	 */
    double  EpochRAAN;			/* RAAN at epoch		 */
    double  epochMeanMotion;		/* Revolutions/day		 */
    double  OrbitalDecay;		/* Revolutions/day^2		 */
    double  EpochArgPerigee;		/* argument of perigee at epoch  */
    double  Eccentricity;
    double  Inclination;
    char    SatName[100];
    int     ElementSet;
    double  BeaconFreq;			/* Mhz, used for doppler calc	 */
    double  MaxPhase;			/* Phase units in 1 orbit	 */
    double  perigeePhase;
    int     NumModes;
    ModeRec Modes[MaxModes];
    int     PrintApogee;
 
    /* Simulation Parameters */
 
    double StartTime,EndTime, StepTime; /* In Days, 1 = New Year	*/
					/*	of reference year	*/
 
    /* Site Parameters */
   
     char SiteName[100];
    double SiteLat,SiteLong,SiteAltitude,SiteMinElev;
 
 
/* List the satellites in kepler.dat, and return the number found */
ListSatellites()
{
    char str[100];
    FILE *InFile;
    char satchar;
    int NumSatellites;

    printf("Available satellites:\n");

    if ((InFile = fopen("kepler.dat","r")) == 0)
	{
	printf(" kepler.dat not found\n");    /* Bitwise Or  "|" */

	exit(-1);
	}

    satchar = 'a';
    NumSatellites = 0;
    while (fgets(str,100,InFile))
 if (sscanf(str,"Satellite: %s",str) == 1)
	    {
	    printf("    %c) %s\n",satchar++,str);
	    NumSatellites++;
	    }

    fclose(InFile);

    return NumSatellites;
}

/* Match and skip over a string in the input file. Exits on failure. */

MatchStr(InFile,FileName,Target)
FILE *InFile;
char *FileName,*Target;
{
    char str[100];

    fgets(str,strlen(Target)+1,InFile);
    if (strcmp(Target,str))
       {
       printf("%s: found \"%s\" while expecting \"%s\n\"",FileName,str,Target);
       exit(-1);
       }
}

GetSatelliteParams()
{
    FILE *InFile;
    char str[100];
    int EpochYear;
   /*    double EpochHour,EpochMinute,EpochSecond;  */
    int found;
    int NumSatellites;
    char satchar;

    NumSatellites = ListSatellites();

    found = 0;

    while (!found)
	{
	printf("Letter or satellite name :");
	gets(SatName);

	if ((InFile = fopen("kepler.dat","r")) == 0)
	    {
	    printf("kepler.dat not found\n");
	    exit(-1);
	    }

	if (strlen(SatName) == 1)
	    {			/* use single character label */
	    satchar = SatName[0];
	    if (satchar < 'a' || satchar > 'a'+NumSatellites-1)
		{
		printf("'%c' is out of range\n",satchar);
		fclose(InFile);
		continue;
		}

	    do	
		do  /* find line beginning with "Satellite: " */
		    fgets(str,100,InFile);
		while (sscanf(str,"Satellite: %s",SatName) != 1);
	     while (satchar-- > 'a');
	    found = 1;
	    }
		
	 else 
	     {
	     while (!found)  /* use satellite name */
		{
		if (! fgets(str,100,InFile))
		    break;	/* EOF */

		if (sscanf(str,"Satellite: %s",str) == 1
		    && strcmp(SatName,str) == 0)
			found = 1;
		}

	    if (!found)
		{
		printf("Satellite %s not found\n",SatName);
		fclose(InFile);
		}
	    }
	}

    BeaconFreq = 146.0;  /* Default value */

    fgets(str,100,InFile);	/* Skip line */

    MatchStr(InFile,"kepler.dat","Epoch time:");
    fgets(str,100,InFile);
    sscanf(str,"%lf",&EpochDay);

    EpochYear = EpochDay / 1000.0;
    EpochDay -= EpochYear*1000.0;
    EpochDay = GetDay(EpochYear,1,EpochDay);

   
    MatchStr(InFile,"kepler.dat","Element set:");
    fgets(str,100,InFile);
    sscanf(str,"%d",&ElementSet);

    MatchStr(InFile,"kepler.dat","Inclination:");
    fgets(str,100,InFile);
    sscanf(str,"%lf",&Inclination);
    Inclination *= RadiansPerDegree;

    MatchStr(InFile,"kepler.dat","RA of node:");
    fgets(str,100,InFile);
    sscanf(str,"%lf",&EpochRAAN);
    EpochRAAN *= RadiansPerDegree;

    MatchStr(InFile,"kepler.dat","Eccentricity:");
    fgets(str,100,InFile);
    sscanf(str,"%lf",&Eccentricity);

    MatchStr(InFile,"kepler.dat","Arg of perigee:");
    fgets(str,100,InFile);
    sscanf(str,"%lf",&EpochArgPerigee);
    EpochArgPerigee *= RadiansPerDegree;

    MatchStr(InFile,"kepler.dat","Mean anomaly:");
    fgets(str,100,InFile);
    sscanf(str,"%lf",&EpochMeanAnomaly);
    EpochMeanAnomaly *= RadiansPerDegree;

    MatchStr(InFile,"kepler.dat","Mean motion:");
    fgets(str,100,InFile);
    sscanf(str,"%lf",&epochMeanMotion);

    MatchStr(InFile,"kepler.dat","Decay rate:");
    fgets(str,100,InFile);
    sscanf(str,"%lf",&OrbitalDecay);

    MatchStr(InFile,"kepler.dat","Epoch rev:");
    fgets(str,100,InFile);
    sscanf(str,"%ld",&EpochOrbitNum);

	while (1)
	    {
	    if (! fgets(str,100,InFile))
		break;	/* EOF */
	    if (strlen(str) <= 2)
		break;	/* Blank line */
	    if (sscanf(str,"Beacon: %lf",&BeaconFreq) == 1)
		break;	/* found it */
	    }


    PrintApogee = (Eccentricity >= 0.3);

    perigeePhase = 0; MaxPhase = 360; /* Default values */
    NumModes = 0;

    if ((InFile = fopen("mode.dat","r")) == 0)
	return;

    found = 0;
    while (!found)
	{
	if (! fgets(str,100,InFile))
	    break;	/* EOF */
	if (sscanf(str,"Satellite: %s",str) == 1
	    && strcmp(SatName,str) == 0)
		found = 1;
	}
	
    if (found)
	{
	MatchStr(InFile,"mode.dat","Perigee phase:");
	fgets(str,100,InFile);
	sscanf(str,"%lf",&perigeePhase);

	MatchStr(InFile,"mode.dat","Max phase:");
	fgets(str,100,InFile);
	sscanf(str,"%lf",&MaxPhase);


	while (fgets(str,100,InFile))
	    {
	    if (sscanf(str,"Mode: %20s from %d to %d",Modes[NumModes].ModeStr,
	    &Modes[NumModes].MinPhase,&Modes[NumModes].MaxPhase) == 3
	      && NumModes < MaxModes)
		  NumModes++;
	    }
	fclose(InFile);
	}



}

 
GetSiteParams()
{
    FILE *InFile;
    char str[100];
 
/*    printf("Site name :");   */
/*    gets(str);    */
/*    strcat(str,".sit");   */
 
    if ((InFile = fopen("Site.dat","r")) == 0)
	{
	printf("%s not found\n",str);
	exit(-1);
	}
 
    fgets(SiteName,100,InFile);
 
    fgets(str,100,InFile);
    sscanf(str,"%lf",&SiteLat);
    SiteLat *= RadiansPerDegree;
 
    fgets(str,100,InFile);
    sscanf(str,"%lf",&SiteLong);
    SiteLong *= RadiansPerDegree;
 
    fgets(str,100,InFile);
    sscanf(str,"%lf",&SiteAltitude);
    SiteAltitude /= 1000;   /* convert to km */
 
    fgets(str,100,InFile);
    sscanf(str,"%lf",&SiteMinElev);
    SiteMinElev *= RadiansPerDegree;
}
 
GetSimulationParams()
{
    double hour,duration;
    int Month,Day,Year;
 
    printf("Start date (UTC) (Month Day Year) :");
    scanf("%d%d%d",&Month,&Day,&Year);
 
    StartTime = GetDay(Year,Month,(double) Day);
 
    printf("Starting Hour (UTC) :");
    scanf("%lf",&hour);
 
 /*  hour = hour + 5;    make up difference for UTC time 
    if ( hour > 24 )  
     hour = hour - 24;     can't have 25 hours corection for over 24 
 */
   StartTime += hour/24;
 
    printf("Duration (Days) :");
    scanf("%lf",&duration);
    EndTime = StartTime + duration;
 
    printf("Time Step (Minutes) :");
    scanf("%lf",&StepTime);
    StepTime /= MinutesPerDay;
}
 
PrintMode(OutFile,Phase)
FILE *OutFile;
{
    int CurMode;
 
    for (CurMode = 0; CurMode < NumModes; CurMode++)
	if ((Phase >= Modes[CurMode].MinPhase
		&& Phase <= Modes[CurMode].MaxPhase)
	      || ((Modes[CurMode].MinPhase > Modes[CurMode].MaxPhase)
		  && (Phase >= Modes[CurMode].MinPhase
			|| Phase <= Modes[CurMode].MaxPhase)))
	    {
	    fprintf(OutFile,"    %s",Modes[CurMode].ModeStr);
	    break;
	    }
}
 
 
main()
{
    double ReferenceOrbit;	/* Floating point orbit # at epoch */
    double CurrentTime; 	/* In Days			   */
    double CurrentOrbit;
    double AverageMotion,CurrentMotion;
    double MeanAnomaly,TrueAnomaly;
    double SemiMajorAxis;
    double Radius;		/* From geocenter		   */
    double RAANPrecession,PerigeePrecession;
    double SSPLat,SSPLong;
    long   OrbitNum,PrevOrbitNum;
    int    Day,PrevDay;
    double Azimuth,Elevation,Range;
    double PrevRange,Velocity,Doppler;
    int    Phase;
    char   FileName[100];
    FILE   *OutFile;
    int    DidApogee;
    double TmpTime,PrevTime;
    int    PrevVisible;
  
    printf("\n");
    printf("%s\n",Str1);
    printf("%s\n",Str2);
    printf("\n\n");

    GetSatelliteParams();
    GetSiteParams();
    GetSimulationParams();
 
    printf("Output file (Filename or * for CRT) :");
    gets(FileName);	/* Skip previous RETURN */
    gets(FileName);
 
 
    if (strlen(FileName) > 0)
	{
	if ((OutFile = fopen(FileName,"w")) == 0)
	    {
	    printf("Can't write to %s\n",FileName);
	    exit(-1);
	    }
	}
      else OutFile = stdout;
 
    fprintf(OutFile,"%s Element Set %d\n",SatName,ElementSet);

    fprintf(OutFile,"%s\n",SiteName);
 
    fprintf(OutFile,"Doppler calculated for freq = %lf MHz\n",BeaconFreq);
 
    SemiMajorAxis = 331.25 * exp(2*log(MinutesPerDay/epochMeanMotion)/3);
    GetPrecession(SemiMajorAxis,Eccentricity,Inclination,&RAANPrecession,
			&PerigeePrecession);
 
    ReferenceOrbit = EpochMeanAnomaly/PI2 + EpochOrbitNum;
 
    PrevDay = -10000; PrevOrbitNum = -10000;
    PrevTime = StartTime-2*StepTime;
 
    BeaconFreq *= 1E6;		/* Convert to Hz */
 
    DidApogee = 0;
 
    /* Start the loop one step early, to init OldRange */
    for (CurrentTime = StartTime-StepTime; CurrentTime <= EndTime;
		CurrentTime += StepTime)
	{
 
	/* Fix these later to correct for drag	*/
	AverageMotion = epochMeanMotion;
	CurrentMotion = epochMeanMotion;
	SemiMajorAxis = 331.25 * exp(2*log(MinutesPerDay/CurrentMotion)/3);
 
	CurrentOrbit = ReferenceOrbit +
			(CurrentTime-EpochDay)*AverageMotion;
	OrbitNum = CurrentOrbit;
 
	MeanAnomaly = (CurrentOrbit-OrbitNum)*PI2;
 
	TmpTime = CurrentTime;
	if (MeanAnomaly < PI)
	    DidApogee = 0;
	if (PrintApogee && !DidApogee && MeanAnomaly > PI)
	    {			/* Calculate Apogee */
	    TmpTime -= StepTime;   /* So we pick up later where we left off */
	    MeanAnomaly = PI;
	    CurrentTime=EpochDay+(OrbitNum-ReferenceOrbit+0.5)/AverageMotion;
	    }
 
	TrueAnomaly = Kepler(MeanAnomaly,Eccentricity);
	Radius = SemiMajorAxis*(1-SQR(Eccentricity))
			/ (1+Eccentricity*cos(TrueAnomaly));
 
	GetSubSatPoint(EpochDay,EpochRAAN,EpochArgPerigee,
	   Inclination,Eccentricity,RAANPrecession,PerigeePrecession,
	   CurrentTime,TrueAnomaly,&SSPLat,&SSPLong);
	GetBearings(SiteLat,SiteLong,SiteAltitude,SSPLat,SSPLong,Radius,
	    &Azimuth,&Elevation,&Range);
 
	Velocity = (Range-PrevRange)/((CurrentTime-PrevTime)*SecondsPerDay);
	PrevRange = Range;
 
	if (Elevation >= SiteMinElev && CurrentTime >= StartTime)
	    {
	    Day = CurrentTime + HalfSecond;
	    if (((double) Day) > CurrentTime+HalfSecond)
		Day -= 1;    /* Correct for truncation of negative values */

	    if (OrbitNum == PrevOrbitNum && Day == PrevDay && !PrevVisible)
		fprintf(OutFile,"\n");	/* Dipped out of sight; print blank */

	    if (OrbitNum != PrevOrbitNum || Day != PrevDay)
		{			/* Print Header */
                fprintf(OutFile,"\n");
                PrintDate(OutFile,Day);
               	fprintf(OutFile,"  ----Orbit # %ld-----\n",OrbitNum);
          	fprintf(OutFile," U.T.C.   Az   El  Doppler  Range");
		fprintf(OutFile,"  Height  Lat  Long  Phase(%3.0lf)\n",
				MaxPhase);
		}
	    PrevOrbitNum = OrbitNum; PrevDay = Day;
            PrintTime(OutFile,CurrentTime + HalfSecond);
 
	    Doppler = -BeaconFreq*Velocity/C;
 
	    fprintf(OutFile,"  %3.0lf  %3.0lf  %6.0lf  %6.0lf",
		Azimuth*DegreesPerRadian,
		Elevation*DegreesPerRadian,Doppler,Range);
 
	    fprintf(OutFile,"  %6.0lf  %3.0lf  %4.0lf",
		Radius - EarthRadius,SSPLat*DegreesPerRadian,
		SSPLong*DegreesPerRadian);
 
	    Phase = (MeanAnomaly/PI2*MaxPhase + perigeePhase);
	    while (Phase < 0)
		Phase += MaxPhase;
	    while (Phase >= MaxPhase)
		Phase -= MaxPhase;
 
	    fprintf(OutFile,"  %4d", Phase);
	    PrintMode(OutFile,Phase);
	    if (PrintApogee && (MeanAnomaly == PI))
		fprintf(OutFile,"    Apogee");
	    fprintf(OutFile,"\n");
	    PrevVisible = 1;
	    }
	 else
	    PrevVisible = 0;	
	if (PrintApogee && (MeanAnomaly == PI))
	    DidApogee = 1;
	PrevTime = CurrentTime;
	CurrentTime = TmpTime;
	}
    fclose(OutFile);
}
   END;