/*
 *  CalculateDrag.c
 *  
 *
 *  Created by Davies, Misty D. (ARC-TI) on 5/10/10.
 *
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <float.h>
#include "dart.h"
#define PI 3.14159

enum varTypes {vInt, vFlt};
enum varTypes vType;

typedef struct Dat {
	int type; // These are the types for the variables
	union aNums { // These are the actual values
		int *inum;
		double *fnum;
	} a;
}dat;

typedef struct Data {
	int numVars; // This is the number of variables
	int numRuns; // This is the number of times we'll use the variables
	char **names; // This is the variable names
	struct Dat *values; // These are the actual values for the variables.  Note that we are assuming numeric values.
}data;

int ParseInputs(struct Data *Input);

// The next two prototypes change for every system
int System(double Pt, double Ps, double Altitude);
double function(double Ma, double Cf, double Cfb, double fb, double CfbTerm, double CfTerm, double L, double d);

double IterSolvePitotMa(double PRat, double gamma);  // A special function just for this one example
double RayleighPitot(double Ma, double gamma); // Ditto.

int main(void)
{
	
	struct Data Inputs;
	int i, Result;
	
	printf("About to parse the inputs.\n");
	ParseInputs(&Inputs);
	
	// There's almost certainly a more elegant solution for this, but as a quick and dirty option this next section is going to be 
	// hardcoded for each System
	double Pt, Ps, Alt;
	FILE *values;
	values=fopen("values.txt","w");
	fprintf(values,"Ma Cf Cfb fb CfbTerm CfTerm L d \n");
	fclose(values);
	for (i=0; i<Inputs.numRuns; i++) {
		Pt=Inputs.values[0].a.fnum[i];
		Ps=Inputs.values[1].a.fnum[i];
		Alt=Inputs.values[2].a.fnum[i];
		
		//print this for now, if it works, we'll actually call the system.
		printf("Inputs: %lf    %lf     %lf\n",Pt, Ps, Alt);
		Result=System(Pt,Ps,Alt);
		printf("Output: %d\n",Result);
	}
	
	// Do all of the frees here.	
	
	
	return 0;
}

int ParseInputs(struct Data *Input) {
	
	FILE *fp;
	char teststring[40];
	int i,j;
	
	// Open the file for reading
	if ((fp = fopen("test_vectors.txt", "r")) == NULL) {
		printf("Can't open test_vectors.txt\n");
		exit(1);
	}
	printf("Opened file\n");
	
	//The first line of the file should be the number of runs
	fscanf(fp,"%s",teststring);
	printf("%s\n",teststring);
	if (strncmp(teststring,"NUM_RUNS:",8)==0) { //strncmp returns 0 if two strings match
		// Read the integer and store it
		fscanf(fp,"%d",&(Input->numRuns));
		printf("Number of runs is %d\n",Input->numRuns);
	}
	else {
		// Something goofy is going on with the file
		printf("Can't parse test_vectors.txt\n");
		exit(1);
	}
	
	// The second line should be the number of variables
	fscanf(fp,"%s",teststring);
	printf("%s\n",teststring);
	if (strncmp(teststring,"NUM_VARS:",8)==0) { //strncmp returns 0 if two strings match
		// Read the integer and store it
		fscanf(fp,"%d",&(Input->numVars));
		printf("Number of vars is %d\n:",Input->numVars);
	}
	else {
		// Something goofy is going on with the file
		printf("Can't parse test_vectors.txt\n");
		exit(1);
	}
	
	// Now we can start to size some variables	
	Input->names=(char **)malloc(Input->numVars*sizeof(char *));
	for (i=0;i<Input->numVars;i++){                
		Input->names[i]=(char *)malloc(80*sizeof(char));
        }
	Input->values=(dat *)malloc(Input->numVars*sizeof(struct Dat));
	
	// First grab the variable names
	fscanf(fp,"%s",teststring);
	printf("%s\n",teststring);
	if (strncmp(teststring,"VARS:",4)==0) { //strncmp returns 0 if two strings match
		// Read the names and store them
		for (i=0;i<Input->numVars;i++) {
			fscanf(fp,"%s",Input->names[i]);
			printf("%s\n",Input->names[i]);
		}
	}
	else {
		// Something goofy is going on with the file
		printf("Can't parse test_vectors.txt\n");
		exit(1);
	}
	
	// Now grab the types
	fscanf(fp,"%s",teststring);
	if (strncmp(teststring,"TYPES:",4)==0) { //strncmp returns 0 if two strings match
		// Read the names and store them
		for (i=0;i<Input->numVars;i++) {
			fscanf(fp,"%s",teststring);
			printf("%s\n",teststring);
			if (strncmp(teststring,"f",1)==0 || strncmp(teststring,"p",1)==0) {
				Input->values[i].type=vFlt;
				Input->values[i].a.fnum=(double *) malloc(Input->numRuns*sizeof(double));
			}
			else if(strncmp(teststring,"i",1)==0) {
				Input->values[i].type=vInt;
				Input->values[i].a.inum=(int *) malloc(Input->numRuns*sizeof(int));
			}
			else {
				// Something goofy is going on with the file
				printf("Can't parse test_vectors.txt\n");
				exit(1);
			}
		}
	}
	else {
		// Something goofy is going on with the file
		printf("Can't parse test_vectors.txt\n");
		exit(1);
	}
	
	
	// Read in the data until we are done
	fscanf(fp,"%s",teststring);
	printf("%s\n",teststring);
	if (strncmp(teststring,"DATA:",4)==0) { //strncmp returns 0 if two strings match
		// Read the names and store them
		printf("matched \n");
		for (i=0;i<Input->numRuns;i++) {
			printf("i is %d\n",i);
			for (j=0; j<Input->numVars;j++) {
				printf("j is %d\n", j);
				printf("The type is %d.\n",Input->values[j].type);
				if (Input->values[j].type==vInt) {
					printf("Putting in an int.\n");
					if (fscanf(fp,"%d",&(Input->values[j].a.inum[i]) < 1)) {
						// Something goofy is going on with the file
						printf("Can't parse test_vectors.txt\n");
						exit(1);
					}
					else
						printf("Should have been a successful read.");	
					printf(" %d\n",Input->values[j].a.inum[i]);
				}
				else {
					printf("Putting in a double.\n");
					if (fscanf(fp,"%lf",&(Input->values[j].a.fnum[i])) < 1) {
						// Something goofy is going on with the file
						printf("Can't parse test_vectors.txt\n");
						exit(1);
					}
					else
						printf("Should have been a successful read.");	
					printf(" %lf\n",Input->values[j].a.fnum[i]);
				}
			}
		}
	}
	else {
		// Something goofy is going on with the file
		printf("Can't parse test_vectors.txt\n");
		exit(1);
	}
	
	// Be nice and close the file
	if (fclose(fp) != 0)
		printf("Error in closing test_vectors.txt\n");

	printf("Done with parsing input file.\n");

	return(0);	
}

int System(double Pt, double Ps, double Altitude)
{
	char replay[50] = "replay";
	char str[10];
	FILE *values;
	double PM1=0;
	double c1,c2;  // These are constants for the kinematic viscosity equation
	double a; // This is the speed of sound
	double Ma; // This is the Mach number
	double Rn;  // This is the Reynolds number
        double K=0.00025;  // This is a surface coefficient, assume smooth matte paint
	double L=60.0;  // This will be the length of our rocket in inches
	double d=12;  // this will be the maximum diameter of the fuselage in inches
	double fb=0;  // Friction on the body
	double Nu=0;  //This is the kinematic viscosity
	double Rnb=0; // This is the compressible Reynolds number;
	double Cf=0, Cfb=0;  // This is the total friction coefficient
	double CfTerm=0, CfbTerm=0; // This is the terminal friction coefficient
	double PRat=Pt/Ps;
        double gamma=1.4;
        double Result;
	double expon;
	
	// First we need to determine what the pitot tube would measure at Mach 1
	PM1=1.893*Ps;
	
	printf("Pt is %lf, Ps is %lf, PM1 is %lf.\n",Pt,Ps,PM1);
	
	// Now we decide whether we are supersonic or subsonic
	if (Pt<PM1)  {
		//This is the subsonic case
                printf("This test case is subsonic.\n");
		expon=(gamma-1)/gamma;
		Ma=sqrt(2/(gamma-1)*(pow(PRat,expon)-1));
		printf("Ma is %lf, the power piece is %lf\n",Ma,pow(PRat,expon)-1);
	}
	else if (Pt>PM1) {
		// This is the supersonic case
		printf("This test case is supersonic.\n");
		Ma=IterSolvePitotMa(PRat,gamma);
	}
	else {
		printf("We are at the critical Mach number.\n");
		Ma=1.0;
	}
	
	
	
	// First we'll need to find the speed of sound in fps, this is a function of altitude
	if (Altitude<=37000) {
		a=-0.004*Altitude+1116.45;
	}
	else if (Altitude>=64000) {
		a=0.0007*Altitude+924.99;
	}
	else {
		a=968.08;
	}
	
	
	// Now we'll need to find the kinematic viscosity in ft^2/s, also a function of altitude
	if (Altitude<=15000) {
		c1=0.00002503;
		c2=0.0;
	}
	else if (Altitude>=30000) {
		c1=0.00004664;
		c2=-0.6882;
	}
	else {
		c1=0.00002760;
		c2=-0.03417;
	}
	Nu=0.000157*exp(c1*Altitude+c2);
	
	// Now the compressible Reynolds number;
	Rn=a*Ma*L/(12*Nu)*(1+0.0283*Ma-0.043*Ma*Ma+0.2107*Ma*Ma*Ma-0.03829*Ma*Ma*Ma*Ma+0.002709*Ma*Ma*Ma*Ma*Ma);
	
	// Now we'll build up all of the parts of the friction coefficient
	// First the incompressible skin friction coefficient
	Cf=0.037036*pow(Rn,-0.155079);
	// Multiply by the compressible skin friction coefficient
	Cf=Cf*(1+0.0283*Ma-0.043*Ma*Ma+0.2107*Ma*Ma*Ma-0.03829*Ma*Ma*Ma*Ma+0.002709*Ma*Ma*Ma*Ma*Ma);
	// Now to find the terminal friction coefficient, start with incompressible skin friction
	CfTerm=1/pow((1.89+1.62*log10(L/K)),2.5);
	// Bring in the compressible skin friction coefficient term
	CfTerm=CfTerm/(1.02044*Ma*Ma);
	
	// First we have to calculate everything for Mach 0.6
	// First the compressible Reynolds number;
	Rnb=a*0.6*L/(12*Nu)*(1+0.0283*0.6-0.043*0.6*0.6+0.2107*0.6*0.6*0.6-0.03829*0.6*0.6*0.6*0.6+0.002709*0.6*0.6*0.6*0.6*0.6);
	
	// Now we'll build up all of the parts of the friction coefficient
	// First the incompressible skin friction coefficient
	Cfb=0.037036*pow(Rnb,-0.155079);
	// Multiply by the compressible skin friction coefficient
	Cfb=Cfb*(1+0.0283*0.6-0.043*0.6*0.6+0.2107*0.6*0.6*0.6-0.03829*0.6*0.6*0.6*0.6+0.002709*0.6*0.6*0.6*0.6*0.6);
	// Now to find the terminal friction coefficient, start with incompressible skin friction
	CfbTerm=1/pow((1.89+1.62*log10(L/K)),2.5);
	// Bring in the compressible skin friction coefficient term
	CfbTerm=CfTerm/(1.02044*0.6*0.6);
	
	printf("About to call the function.\n");
	Result=function(Ma, Cf, Cfb, fb, CfbTerm, CfTerm, L, d);
	
	printf("Run returns %lf\n", Result);
	values=fopen("values.txt","a");
	//fprintf(values,"Ma Cf Cfb fb CfbTerm CfTerm L d \n");
	
	if (values==NULL) {
	}
	else {
		fprintf(values,"%lf %lf %lf %lf %lf %lf %lf %lf\n",Ma, Cf, Cfb, fb, CfbTerm, CfTerm, L, d);
		fclose(values);
	}
	
	dart_generateSummary("summary.txt", CSyntax);
	dart_solve(1);
	return 0;
}

double function(double Ma, double Cf, double Cfb, double fb, double CfbTerm, double CfTerm, double L, double d) {
	
	// This function calculates the drag coefficient, given the current altitude
	// and Mach number.  To make the problem interesting, we'll do the math and unit 
	// conversions here
	
	double Cd=0, Cdb=0; // This is the drag coefficent;
	double db=4;  // The diameter at the tail
	double Sb=4070; // This is the wetted surface area, in inches squared;
	double Ke=0;  // Excrescencies drag increment

	

	
	if (Cf<=CfTerm)
		Cf=CfTerm;
	
	Cd=Cf*(1+60/pow((L/d),3)+0.0025*(L/d))*4*Sb/(PI*d*d);
	
	
	if (Ma<0.78)
		Ke=0.78;
	else if (Ma>1.04)
		Ke=0.0002*Ma*Ma-0.0012*Ma+0.0018;
	else
		Ke=-0.4501*Ma*Ma*Ma*Ma+1.5954*Ma*Ma*Ma-2.1062*Ma*Ma+1.2288*Ma-0.26717;
	
	Cf=Cf+Cd+Ke*4*Sb/(PI*d*d);
	
	Cd=Cd+Ke*4*Sb/(PI*d*d);
	
	// Now add the Base Drag Coefficient
	
	if (Ma<0.6)
		Cd=Cd+0.029*pow((db/d),3)/sqrt(Cf);
	else {

		
		if (Cfb<=CfbTerm)
			Cfb=CfbTerm;
		
		Cdb=Cfb*(1+60/pow((L/d),3)+0.0025*(L/d))*4*Sb/(PI*d*d);
		
		// now switch between transonic, supersonic, and hypersonic
		if (Ma<1.0)
			fb=1.0+215.8*pow((Ma-0.6),6.0);
		else if (Ma>2.0)
			fb=0.297*pow((Ma-2.0),3.)-0.7937*pow((Ma-2.0),2)-0.1115*(Ma-2.0)+1.64006;
		else
			fb=2.0881*pow((Ma-1.),3)-3.7938*pow((Ma-1),2)+1.4618*(Ma-1.0)+1.883917;
		
		Cd=Cd+Cdb*fb;
		
	}
	
	//dart_printReplayData("replay.txt", "params.dat", DartAppendMode); 
	return Cd;
	
}

double IterSolvePitotMa(double PRat, double gamma) {
	// This function finds Ma for supersonic flow.  It requires an iterative solve.
	double tol=0.001;  // We really want to nail down the Mach number
	
	
	
	double MaTry=1.2;  // This is our initial guess for the Mach number
	int converged=0; //  A boolean to see if we've converged
	int timesthru=500; //  An integer max number of times we'll loop
	int thistime=0; // A counter to see how many times we've gone around
	
	double MaNew=0;
	
	// We'll also need the derivative.
	double DT1=(2*gamma)/(1+gamma);
	
	double Fx=0;
	double Fxp=0;

	double h=sqrt(DBL_EPSILON);
	
	// Use Newton's Method to find the solution
	while ((converged==0) && (thistime<=timesthru)) {
		// find the solution at the current value of MaTry
		Fx=RayleighPitot(MaTry,gamma)-PRat;
		Fxp=(RayleighPitot(MaTry+h,gamma)-RayleighPitot(MaTry-h,gamma))/(2*h);
		printf("Fx is %lf, Fxp is %lf\n",Fx,Fxp);
		MaNew= MaTry-Fx/Fxp;
	
		if (MaNew<1.0)
			MaNew=1.0;
	
		printf("MaNew is %lf, MaTry is %lf\n",MaNew,MaTry);
		//See if we converged
		if (fabs(MaNew-MaTry)/(fabs(MaTry)+1)<=tol)
			converged=1;
		
		//Update the counter and MaTry
		thistime=thistime+1;
		MaTry=MaNew;
		
		if (thistime>=timesthru)
			printf("WARNING:  Iterative solver timed out.\n");
		
	}
	
	return MaNew;
	
	
	
}

double RayleighPitot(double Ma, double gamma) {



	// Precalculate the terms that won't change
	double T2Denom=(gamma+1);
	double T1Num=T2Denom*T2Denom;
	double T2Num1=1-gamma;
	double T2Num2=2*gamma;
	double T1Power=gamma/(gamma-1);
	double T1Denom1=4*gamma;
	double T1Denom2=2*(gamma-1);
	double Prat;

	Prat=pow(((T1Num*Ma*Ma)/(T1Denom1*Ma*Ma-T1Denom2)),T1Power)*(T2Num1+T2Num2*Ma*Ma)/T2Denom;

	return Prat;
}