#include "asa_user.h"
#include <math.h>
 #define LaborFactor (1.0)  // labor factor on stuff, may divide into subfactors <reg>//Parameter// set to 1 for experimental 
//Note:  CONSTANT denotes a questionable value used as a define
//       Constant denotes a real constant, grey areas marked examine them later CONSTANT
//Parameter is a parameter to the simulation, perhaps these should come
// from somewhere else, but it might be preferable to just group them
// together in this file for now
//#define CarbonTankDensity 2202.7//Constant
// from http://www.lmco.com/manned/helium.html burst pressure of 7,200psi
// .25in/(16in*Pi) = 0.0049736//Constant
//#define CarbonTankThicknessPerCircumference (0.0049736/4) // bp 1,800psi//Constant
void setupExpArray();
double CalcISP(double, double);
static inline double Throttle(double, double, double, double , double);
int ExpArrayInitialized = 0;
#define MaxSECS 140002
#define DegreesToRadians(x) (x*Pi/180.)//CONSTANT
#define RadiansToDegrees(x) (x*180./Pi)//CONSTANT
#define NewtonsPerLbf 4.448222  //CONSTANT
#define NewtonsToLbf(x)  (x/NewtonsPerLbf) 
#define PascalsPerPsi 6894.75 //Constant
#define PascalsToPsi(x) (x/PascalsPerPsi)
#define PsiToPascals(x) (x*PascalsPerPsi)
#define AluminumDensity (2702.0) //C//CONSTANTonstant
#define PET 0
#define PETDensity       (1370.)////CONSTANTConstant
#define PETStrength       (3.44379e8)////CONSTANTConstant
#define PETPrice	   (5.) // Constant guess
#define FUZZ .0001//CONSTANT
#define Dyneema 1
#define DyneemaDensity 970 //  this IS right http://dns.toyobo.co.jp/e/seihin/dn/dyneema/karui.htm//Constant//CONSTANT
#define DyneemaStrength ((3e9)/4)	// 3 giga pascals corrected to a working strength with a safety factor of 2 thus a total factor of 4  http://dns.toyobo.co.jp/e/seihin/dn/dyneema/kyoudo.htm//Constant//CONSTANT
#define DyneemaPrice (40*2.2)	//$40/lb but this was for aramid, not dyneema//Parameter//Constant
#define DcWater (.97*.62)//CONSTANT
#define GEarth 9.80665//CONSTANT
#define EarthRadius 6377500.//CONSTANT
#define EarthEquatorialRotationalVelocity (2*Pi*EarthRadius/(24* 3600))//CONSTANT
#define InitialAltitude 0.
//#define InitialAltitude 2000.
//#define InitialAltitude 4000.
//#define InitialAltitude 200. // brutus!
//#define Latitude 28.24// 28.24 is cape canaveral  //Parameter
//#define Latitude 00.//  equatorial orbit
#define Latitude 90.//  polar orbit
//#define Latitude 47.//  neuchatel
//#define Latitude 46.//  brutus
#define EarthRotationalVelocity (EarthEquatorialRotationalVelocity *cos(DegreesToRadians(Latitude)))
//#define FuelPress PsiToPascals(5000.)//Parameter
#define ChamberPressConstant (PsiToPascals(10000))//Parameter
#define NominalMaxISP  400. //Parameter
#define Fudge .8//Parameter
//#define Fudge .90//Parameter
#define CharisticISP 214.3 //from sutton for LOX&RP1 with a CF_Thrust_coefficient of 1.4  pp194, 60//Parameter

#define KSpecificHeatRatio (1.24)  //Parameter sutton table 6-4 p 194
// KSpecificHeatRatio is the specific reat ratio as defined by KSpecificHeadRatio
//I won't bother to simplify this, since it is all constants, and
//the compiler can deal with completely!!! how nice for expressive
//purposes!
//#define ExpansionRatio (x[4])
//#define ExpansionRatio 25. //Parameter corresponds to optimal at 15psi & 4000psi with k = 1.3 //Parameter
//#define ExpansionRatio 40. //Parameter
#define AproxISP 330
#define AproxGThrottle fmin (1.,((MaxDesiredG*VehicleMass/AproxISP)/MassFlowAllowed))
#define ChamberPress (ChamberPressConstant * AproxGThrottle)
//#define CalcMassFlow (Throttle(Y,a,b,c,d)*(OutOfFuel?0.:(fmin(MaxDesiredG*VehicleMass/CalcISP(AirPressure,ChamberPress) ,(fmin(MassFlowAllowed,(TankCrossSection/ExpansionRatio)*(ChamberPress/CalcISP(AirPressure,ChamberPress))))))))// stop burning when out of fuel(AirPressure)
//#define CalcMassFlow (Throttle(Y,a,b,c,d)*(OutOfFuel?0.:(fmin(MaxDesiredG*VehicleMass/CalcISP(AirPressure,ChamberPress) ,(fmin(MassFlowAllowed,(TankCrossSection/ExpansionRatio)*(ChamberPress/CalcISP(AirPressure,ChamberPress))))))))// stop burning when out of fuel(AirPressure)
#define CalcMassFlow ((OutOfFuel?0.:(fmin(MaxDesiredG*VehicleMass/CalcISP(AirPressure,ChamberPress) ,(fmin(MassFlowAllowed,(TankCrossSection/ExpansionRatio)*(ChamberPress/CalcISP(AirPressure,ChamberPress))))))))// stop burning when out of fuel(AirPressure)
#define CalcThrust (CalcMassFlow*CalcISP(AirPressure,ChamberPress)*GEarth)

#define ThrustAndMassToAcc ((CalcThrust-Drag)/VehicleMass)
#define FuelVolume (FuelMass/FuelDensity)
#define FuelUnitCost 500.0	// about $500/cubic meter of LOX (or fuel)//Parameter/CONSTANTrameter/
#define FuelCost (FuelUnitCost*TankVolume)
#define FuelHeight  (FuelVolume/TankCrossSection)
#define TankBottomPress (TankPress+Acc*(FuelHeight*FuelDensity))//do we have teh units right here? zzz
//#define TankPress PsiToPascals(100.0)//Parameter

#define VehicleMass  (VehicleDryMass+FuelMass+PayLoad)
#define MeritVehicleCost (FuelCost+LaborFactor*(TankCost+EngineCost))
#define VehicleCost (FuelCost+(TankCost+EngineCost))
#define TankCost TankMass*MaterialPrice
#define EngineCostPerWeight 50.//Parameter
#define EngineCost  (fmax(EngineCostPerWeight*0.,(EngineMass *EngineCostPerWeight)))
#define UninsuredLaunchCost (MeritVehicleCost)
#define RelaunchInsurance (UninsuredLaunchCost*InsuranceRate)
#define PayLoadInsurance (PayLoadValue*InsuranceRate)
#define PayLoadValue Price	// Guess that payload is about same value as launch
//#define LaunchCost (UninsuredLaunchCost+RelaunchInsurance+PayLoadInsurance)
#define LaunchCost (UninsuredLaunchCost)
#define InsuranceRate .4//Parameter
/* *****************************************************
// I use fmin fmax & imin imax inlinw functions here for efficiency.  Most calls
// to these are to expensive functions, do I don't want double evaluation.
********************************************************** */
//#define min(x,y) ((x)<(y)?(x):(y))
//#define imin(x,y) ((x)<(y)?(x):(y)) // for integers
static inline int imin(int x,int y){
  return x<y?x:y;
}
//#define min(x,y) fmin((x),(y)) 
static inline double fmin(double xx,double yy){  // for doubles
	return	xx > yy ? yy : xx;
}
#define max(x,y) ((x)>(y)?(x):(y))
//#define max(x,y) fmax((x),(y))
static inline double fmax(double xx,double yy){
	return	xx > yy ? xx : yy;
}

#define MilesToM(x) (x*1609.344)//CONSTANT
//#define OrbitAltitude             482800.//Parameter  //300 miles LEO
#define OrbitAltitude             MilesToM(110)//Parameter // 110 miles LEO
//#define OrbitAltitude             MilesToM(250000)//Parameter // 110 miles LEO
//#define OrbitAltitude  100000.  // X prize altitude
//#define OrbitAltitude             35786000.  //geosync//Parameter
//#define AltitudeOnly TRUE // only try to match altitude, current design of trajedctory is crude, this helps.
#define AltitudeOnly  FALSE
#define MadeOrbit  ((AltitudeOnly?1:VelocityX>OrbitalVelocityAtAltitude)&&(Y>OrbitAltitude))//
#define Price (MadeOrbit?PayLoad*(2000.0*2.2):0)////Parameter
#define Profit (Price-LaunchCost)
//#define Merit (MadeOrbit?(Profit/UninsuredLaunchCost)-aspectBogosityFactor-GBogosityFactor:(fmin(OrbitAltitude,fmax(1,Y)) *fmin(OrbitalVelocityAtAltitude,fmax(1,VelocityX))-10e10))//
#define Merit (MadeOrbit?(Profit/UninsuredLaunchCost)+((TankPress-10)/MaxFPress<.6 ?(-100.0*MaxFPress/TankPress-1000):0)-aspectBogosityFactor:(fmin(OrbitAltitude,fmax(1,Y)) *fmin(OrbitalVelocityAtAltitude,fmax(1,VelocityX))-10e10))//
//#define Merit (MadeOrbit?(Profit/UninsuredLaunchCost):(fmin(OrbitAltitude,fmax(1,Y)) *fmin(OrbitalVelocityAtAltitude,fmax(1,VelocityX))-10e10))//

//#define Merit (MadeOrbit?10e10-VehicleCost:(fmax(OrbitAltitude,Y)*fmin(OrbitalVelocityAtAltitude,VelocityX))-10e10)//
#define MinimalMaterialThickness .00025//Parameter
#define Pi (3.141592654)//Constant  !! it sure is!!
//	#define TankRadius (TankCircumference/(2.0*Pi))
#define	TankDensity  MaterialDensity
//#define MaxDesiredG 5.//Parameter
#define MinDesiredG 0.//Parameter
#define MaxDesiredAcc (MaxDesiredG*GEarth)
#define MinDesiredAcc (MinDesiredG*GEarth)
#define TankLinerDensity AluminumDensity
#define TankLinerThickness .00025//Parameter
#define	EngineThickness (AvgTankThickness*2.0)// just bogus
#define	EngineLength (4.0)//CONSTANT likewise bogus
#define	EngineDensity (2*AluminumDensity)
#define  ThrustWeightRatioOfEngine (500.) // 100 high normal range we might be 300 or 1000//Parameter
#define MinEngineMass .001 //Parameter
#define NominalEngineMass (MaxThrust/(ThrustWeightRatioOfEngine*GEarth))
// propane c3h8 + 7O2    44 + 160 = 204
#define OxidiserFractionByMass .7843
#define FuelFractionByMass .2157
#define OxidiserDensity 1140//Parameter
#define PropaneDensity 594 // at 200K and >15psi //Parameter
#define PressureTankFraction .3 // the ammount of empty tankage we leave for pressurization, later we will add this in as a function of tank pressure and let the anealer reduce the tank pressure and thickness, leaving some empty space for gass to pressureize the system.  A little space will give a large pressure swing,a lot will give a smaller swing but more tankage, my guess is that the optimum is between .5  and 1-1/e 
#define FuelSpecificVolume  (OxidiserFractionByMass/OxidiserDensity+FuelFractionByMass/PropaneDensity)
//  .630 + .402 = 1.03   1/x = .968
#define FuelDensity  ((1/FuelSpecificVolume))
#define CoefficientOfDrag .20//CONSTANT // set a bit high on purpose for now
#define Drag (CoefficientOfDrag/2.0)*AirDensity*(Velocity)*(Velocity)*TankCrossSection // what units is Drag ? zzz  newtons??
  static double airdens[300];
	
/* ***************************************************************
//
//  so we hackup a function to replace the #define we were using.
//  it interpolates between the 1km resolution table.  This is historical
//  there was a table here, and efficiency is a bit important so we keep it.
//
*****************************************************************/
#define AirDensity  air_(Y)
static inline double air_(double y){ 
  int i;
  double a = 0.;
  
  int km;
	i = y/1000.;
	if(i>200){
		i = 200;
	}
	if (i<0){
		i = 0;
	}	
	km = y/1000.;
	a =  (y - km * 1000.)/1000.;
	  return airdens[i] - a*(airdens[i]  - airdens[i+1]);
}


#define AirPressure PsiToPascals((AirDensity/airdens[0])* 15.)
//#define AirPressure ((AirDensity/airdens[0])* 15.)

extern int _flsbuf (unsigned int, FILE *);

  int NSec = 0;
static double ExpansionRatio;
double
rocket_cost_function (double *x,
	       double *parameter_lower_bound,
	       double *parameter_upper_bound,
               double *cost_tangents,
               double *cost_curvature,
               ALLOC_INT * parameter_dimension,
               int *parameter_int_real,
               int *cost_flag,
               int *exit_code,
               USER_DEFINES * USER_OPTIONS)

{
void checkMR(double,double,double);
  static LONG_INT funevals = 0;

  double TankHeight=x[0];
  double TankPseudoVolume=x[1];
  double MassFlowAllowed=x[2];
  double FlatVelocity=x[3];
  //ExpansionRatio = x[4];  
  double a = x[5];
  double b = x[6];
  double c = x[7];
  double d = x[8];
	double MaxDesiredG = x[9];
	double TankPressE = x[10];
	double TankPress = PsiToPascals(TankPressE);
  int MaterialType = 0;//0 PET 1 Dyneema + liner

double  TankCrossSection = TankPseudoVolume / TankHeight;
double TankRadius = sqrt (TankCrossSection /Pi);
//double TankRadius = TankCircumference/(2.0*Pi);
  double TankCircumference = 2* Pi * TankRadius;
//double TankCrossSection =  ((TankRadius*TankRadius)*Pi);
//double TankVolume =  (TankHeight*TankCrossSection+(4/3)*Pi*(TankRadius*TankRadius*TankRadius)*2);
double TankVolume =  (TankHeight*TankCrossSection+(4/3)*Pi*(TankRadius*TankRadius*TankRadius));// common bulkhead
/* the following Tank properties depend on material selection so they require initialization after that is picked */
 double TankMass = 3e15;
  double TankStructureMass = 4e14;
  double TankLinerMass = 0;
  double EngineMass = 0;
  double VehicleDryMass = 0;
	double sinOrbit;
	double yComponent;
	double xComponent;

/*   */



double PayLoad = 100. ;  // //Parameter
//double PayLoad = x[12] ;  // //Parameter
//double PayLoad = 1. ;  // //Parameter
 double MaxThrust=NominalMaxISP*MassFlowAllowed*GEarth;
// Initial conditions
//  double TankHeight = TankCircumference if TankHeight<TankCircumference;
  double Acc = GEarth;
  double VelocityX=EarthRotationalVelocity,VelocityY=0,Velocity=0,FPress=0;
  double Y = InitialAltitude;
  double X = 0;
  double MFR = 0;
  double Theta = Pi/2.0;
  double THRUSTSUM=0;
  double DeltaV = 0;
  double VMass,ISP,Thrust=0;
  double MaxAcc=0,MinAcc=10000000.;
  double MinAccVelocity=0,MinAccAltitude=0;
  double MaxFPress=0;
  double MaxDrag=0,MaxDragAlt=0, MaxDragV=0;
  double MaxISPAltitude = 0;
  double MaxISP = 0;
  double MinISP = 999999;
  double MaxY = 0, MaxV = 0;
  double FuelMass = 0;
  double StartFuelMass=0;
  double OrbitalVelocityAtAltitude =  EarthRadius * sqrt(GEarth/(EarthRadius + OrbitAltitude));
static double MaxMerit=-10e18;
static int testNumber = 0;
static int printedTestNumber = 0;	
  double GBogosityFactor = 0;
  double TakeOffThrust = 0;
 double FractionOfDesiredHeightCalculated = 0;
 double MaterialPrice = 0;
 double MaterialStrength = 0;
 double MaterialDensity = 0;
 char * MaterialName = NULL;
 int RequiresLiner = 0;
	int OutOfFuel = 0;
	int Burnout = 0;
	double BurnoutAltitude = 0;
	double BurnoutVelocity = 0;
 static int debugPrint = 0;
  double TotalDrag = 0; 
double aspectBogosityFactor;
double	ActualTankThickness  = 3.91e-17;
double AlternateTankThickness = 3.7e-22;
double Xarray[MaxSECS];
double Yarray[MaxSECS];
double YVelocityarray[MaxSECS];
double XVelocityarray[MaxSECS];
double ISParray[MaxSECS];
double Phiarray[MaxSECS];
double Thrustarray[MaxSECS];

int i;
if(airdens[0]!=1.2250){
/*
int i;
     for(i=0;i<300;i++){
	airdens[i]=1.2250*exp(i*(-1.42e-1));
	//printf("airdens[%d] = %f;\n",i, 1.2250*exp(i*(-1.42e-1)));
	//printf("air[%d] = %f;\n",i, exp(i*(-1.42e-1)));
     }
*/


airdens[0] = 1.225;
airdens[1] = 1.06284;
airdens[2] = 0.922139;
airdens[3] = 0.800068;
airdens[4] = 0.694156;
airdens[5] = 0.602264;
airdens[6] = 0.522537;
airdens[7] = 0.453364;
airdens[8] = 0.393349;
airdens[9] = 0.341278;
airdens[10] = 0.2961;
airdens[11] = 0.256902;
airdens[12] = 0.222894;
airdens[13] = 0.193388;
airdens[14] = 0.167787;
airdens[15] = 0.145576;
airdens[16] = 0.126305;
airdens[17] = 0.109585;
airdens[18] = 0.0950779;
airdens[19] = 0.0824916;
airdens[20] = 0.0715714;
airdens[21] = 0.0620969;
airdens[22] = 0.0538766;
airdens[23] = 0.0467445;
airdens[24] = 0.0405565;
airdens[25] = 0.0351877;
airdens[26] = 0.0305296;
airdens[27] = 0.0264881;
airdens[28] = 0.0229817;
airdens[29] = 0.0199394;
airdens[30] = 0.0172998;
airdens[31] = 0.0150097;
airdens[32] = 0.0130227;
airdens[33] = 0.0112988;
airdens[34] = 0.00980308;
airdens[35] = 0.00850536;
airdens[36] = 0.00737943;
airdens[37] = 0.00640255;
airdens[38] = 0.00555499;
airdens[39] = 0.00481962;
airdens[40] = 0.00418161;
airdens[41] = 0.00362805;
airdens[42] = 0.00314778;
airdens[43] = 0.00273108;
airdens[44] = 0.00236954;
airdens[45] = 0.00205586;
airdens[46] = 0.00178371;
airdens[47] = 0.00154759;
airdens[48] = 0.00134272;
airdens[49] = 0.00116497;
airdens[50] = 0.00101075;
airdens[51] = 0.000876951;
airdens[52] = 0.000760862;
airdens[53] = 0.00066014;
airdens[54] = 0.000572751;
airdens[55] = 0.000496931;
airdens[56] = 0.000431148;
airdens[57] = 0.000374073;
airdens[58] = 0.000324554;
airdens[59] = 0.00028159;
airdens[60] = 0.000244313;
airdens[61] = 0.000211971;
airdens[62] = 0.000183911;
airdens[63] = 0.000159565;
airdens[64] = 0.000138442;
airdens[65] = 0.000120115;
airdens[66] = 0.000104215;
airdens[67] = 9.04187e-05;
airdens[68] = 7.84492e-05;
airdens[69] = 6.80642e-05;
airdens[70] = 5.90539e-05;
airdens[71] = 5.12365e-05;
airdens[72] = 4.44538e-05;
airdens[73] = 3.85691e-05;
airdens[74] = 3.34634e-05;
airdens[75] = 2.90335e-05;
airdens[76] = 2.51901e-05;
airdens[77] = 2.18555e-05;
airdens[78] = 1.89623e-05;
airdens[79] = 1.64521e-05;
airdens[80] = 1.42742e-05;
airdens[81] = 1.23846e-05;
airdens[82] = 1.07451e-05;
airdens[83] = 9.32269e-06;
airdens[84] = 8.08857e-06;
airdens[85] = 7.01781e-06;
airdens[86] = 6.0888e-06;
airdens[87] = 5.28277e-06;
airdens[88] = 4.58345e-06;
airdens[89] = 3.9767e-06;
airdens[90] = 3.45027e-06;
airdens[91] = 2.99352e-06;
airdens[92] = 2.59725e-06;
airdens[93] = 2.25343e-06;
airdens[94] = 1.95512e-06;
airdens[95] = 1.6963e-06;
airdens[96] = 1.47175e-06;
airdens[97] = 1.27692e-06;
airdens[98] = 1.10788e-06;
airdens[99] = 9.61223e-07;
airdens[100] = 8.33978e-07;
airdens[101] = 7.23577e-07;
airdens[102] = 6.27791e-07;
airdens[103] = 5.44684e-07;
airdens[104] = 4.7258e-07;
airdens[105] = 4.1002e-07;
airdens[106] = 3.55742e-07;
airdens[107] = 3.0865e-07;
airdens[108] = 2.67791e-07;
airdens[109] = 2.32341e-07;
airdens[110] = 2.01584e-07;
airdens[111] = 1.74899e-07;
airdens[112] = 1.51746e-07;
airdens[113] = 1.31658e-07;
airdens[114] = 1.14229e-07;
airdens[115] = 9.91077e-08;
airdens[116] = 8.59879e-08;
airdens[117] = 7.46049e-08;
airdens[118] = 6.47288e-08;
airdens[119] = 5.61601e-08;
airdens[120] = 4.87257e-08;
airdens[121] = 4.22755e-08;
airdens[122] = 3.66791e-08;
airdens[123] = 3.18236e-08;
airdens[124] = 2.76108e-08;
airdens[125] = 2.39557e-08;
airdens[126] = 2.07845e-08;
airdens[127] = 1.80331e-08;
airdens[128] = 1.56459e-08;
airdens[129] = 1.35747e-08;
airdens[130] = 1.17777e-08;
airdens[131] = 1.02186e-08;
airdens[132] = 8.86585e-09;
airdens[133] = 7.6922e-09;
airdens[134] = 6.67392e-09;
airdens[135] = 5.79043e-09;
airdens[136] = 5.0239e-09;
airdens[137] = 4.35884e-09;
airdens[138] = 3.78182e-09;
airdens[139] = 3.28119e-09;
airdens[140] = 2.84683e-09;
airdens[141] = 2.46997e-09;
airdens[142] = 2.143e-09;
airdens[143] = 1.85931e-09;
airdens[144] = 1.61318e-09;
airdens[145] = 1.39963e-09;
airdens[146] = 1.21435e-09;
airdens[147] = 1.05359e-09;
airdens[148] = 9.1412e-10;
airdens[149] = 7.9311e-10;
airdens[150] = 6.88119e-10;
airdens[151] = 5.97027e-10;
airdens[152] = 5.17993e-10;
airdens[153] = 4.49422e-10;
airdens[154] = 3.89928e-10;
airdens[155] = 3.3831e-10;
airdens[156] = 2.93525e-10;
airdens[157] = 2.54668e-10;
airdens[158] = 2.20956e-10;
airdens[159] = 1.91706e-10;
airdens[160] = 1.66328e-10;
airdens[161] = 1.4431e-10;
airdens[162] = 1.25206e-10;
airdens[163] = 1.08632e-10;
airdens[164] = 9.4251e-11;
airdens[165] = 8.17742e-11;
airdens[166] = 7.0949e-11;
airdens[167] = 6.15569e-11;
airdens[168] = 5.34081e-11;
airdens[169] = 4.6338e-11;
airdens[170] = 4.02038e-11;
airdens[171] = 3.48817e-11;
airdens[172] = 3.02641e-11;
airdens[173] = 2.62578e-11;
airdens[174] = 2.27818e-11;
airdens[175] = 1.9766e-11;
airdens[176] = 1.71494e-11;
airdens[177] = 1.48792e-11;
airdens[178] = 1.29095e-11;
airdens[179] = 1.12005e-11;
airdens[180] = 9.71783e-12;
airdens[181] = 8.43139e-12;
airdens[182] = 7.31526e-12;
airdens[183] = 6.34687e-12;
airdens[184] = 5.50668e-12;
airdens[185] = 4.77771e-12;
airdens[186] = 4.14525e-12;
airdens[187] = 3.5965e-12;
airdens[188] = 3.1204e-12;
airdens[189] = 2.70733e-12;
airdens[190] = 2.34893e-12;
airdens[191] = 2.03799e-12;
airdens[192] = 1.7682e-12;
airdens[193] = 1.53413e-12;
airdens[194] = 1.33104e-12;
airdens[195] = 1.15484e-12;
airdens[196] = 1.00196e-12;
airdens[197] = 8.69325e-13;
airdens[198] = 7.54245e-13;
airdens[199] = 6.54399e-13;
airdens[200] = 5.6777e-13;
airdens[201] = 4.9261e-13;
airdens[202] = 4.27399e-13;
airdens[203] = 3.7082e-13;
airdens[204] = 3.21731e-13;
airdens[205] = 2.79141e-13;
airdens[206] = 2.42189e-13;
airdens[207] = 2.10128e-13;
airdens[208] = 1.82312e-13;
airdens[209] = 1.58177e-13;
airdens[210] = 1.37238e-13;
airdens[211] = 1.19071e-13;
airdens[212] = 1.03308e-13;
airdens[213] = 8.96324e-14;
airdens[214] = 7.7767e-14;
airdens[215] = 6.74723e-14;
airdens[216] = 5.85404e-14;
airdens[217] = 5.07909e-14;
airdens[218] = 4.40673e-14;
airdens[219] = 3.82337e-14;
airdens[220] = 3.31724e-14;
airdens[221] = 2.87811e-14;
airdens[222] = 2.49711e-14;
airdens[223] = 2.16654e-14;
airdens[224] = 1.87974e-14;
airdens[225] = 1.6309e-14;
airdens[226] = 1.415e-14;
airdens[227] = 1.22769e-14;
airdens[228] = 1.06517e-14;
airdens[229] = 9.24162e-15;
airdens[230] = 8.01823e-15;
airdens[231] = 6.95678e-15;
airdens[232] = 6.03585e-15;
airdens[233] = 5.23683e-15;
airdens[234] = 4.54359e-15;
airdens[235] = 3.94211e-15;
airdens[236] = 3.42026e-15;
airdens[237] = 2.96749e-15;
airdens[238] = 2.57466e-15;
airdens[239] = 2.23383e-15;
airdens[240] = 1.93812e-15;
airdens[241] = 1.68155e-15;
airdens[242] = 1.45895e-15;
airdens[243] = 1.26582e-15;
airdens[244] = 1.09825e-15;
airdens[245] = 9.52864e-16;
airdens[246] = 8.26725e-16;
airdens[247] = 7.17284e-16;
airdens[248] = 6.22331e-16;
airdens[249] = 5.39948e-16;
airdens[250] = 4.6847e-16;
airdens[251] = 4.06455e-16;
airdens[252] = 3.52649e-16;
airdens[253] = 3.05966e-16;
airdens[254] = 2.65462e-16;
airdens[255] = 2.30321e-16;
airdens[256] = 1.99831e-16;
airdens[257] = 1.73378e-16;
airdens[258] = 1.50426e-16;
airdens[259] = 1.30513e-16;
airdens[260] = 1.13236e-16;
airdens[261] = 9.82458e-17;
airdens[262] = 8.52401e-17;
airdens[263] = 7.39562e-17;
airdens[264] = 6.41659e-17;
airdens[265] = 5.56717e-17;
airdens[266] = 4.8302e-17;
airdens[267] = 4.19078e-17;
airdens[268] = 3.63601e-17;
airdens[269] = 3.15468e-17;
airdens[270] = 2.73707e-17;
airdens[271] = 2.37474e-17;
airdens[272] = 2.06037e-17;
airdens[273] = 1.78762e-17;
airdens[274] = 1.55098e-17;
airdens[275] = 1.34566e-17;
airdens[276] = 1.16753e-17;
airdens[277] = 1.01297e-17;
airdens[278] = 8.78875e-18;
airdens[279] = 7.62531e-18;
airdens[280] = 6.61588e-18;
airdens[281] = 5.74008e-18;
airdens[282] = 4.98021e-18;
airdens[283] = 4.32094e-18;
airdens[284] = 3.74894e-18;
airdens[285] = 3.25266e-18;
airdens[286] = 2.82208e-18;
airdens[287] = 2.44849e-18;
airdens[288] = 2.12436e-18;
airdens[289] = 1.84314e-18;
airdens[290] = 1.59915e-18;
airdens[291] = 1.38746e-18;
airdens[292] = 1.20379e-18;
airdens[293] = 1.04443e-18;
airdens[294] = 9.06171e-19;
airdens[295] = 7.86213e-19;
airdens[296] = 6.82135e-19;
airdens[297] = 5.91835e-19;
airdens[298] = 5.13489e-19;
airdens[299] = 4.45514e-19;
}

	if(!ExpArrayInitialized){
		ExpArrayInitialized = 1;
		setupExpArray();
	}
   NSec = 0;
  ExpansionRatio = x[4];  
	switch (MaterialType){
              case PET: //PET
                  MaterialPrice =PETPrice;	
                  MaterialStrength = PETStrength;
                  MaterialDensity = PETDensity;
                  MaterialName = "PET";
		  RequiresLiner = FALSE;
                  break;
               case Dyneema://Dyneema
                  MaterialPrice =DyneemaPrice;	
                  MaterialStrength = DyneemaStrength;
                  MaterialDensity = DyneemaDensity;
                  MaterialName = "Dyneema";
		  RequiresLiner = TRUE;
                  break;

}

 AlternateTankThickness =  ((TankPress*2*TankRadius)/MaterialStrength); // standard formula stress = Pressure*Diameter 
	ActualTankThickness = fmax(MinimalMaterialThickness,AlternateTankThickness);
/* initialization stuff that depends on material  */
/* these used to be macros and rere recalculated 1000s of times */
   FuelMass = TankVolume*FuelDensity*(1-PressureTankFraction);
  StartFuelMass=FuelMass;
        EngineMass = (fmax(MinEngineMass,NominalEngineMass));
	TankLinerMass =(TankHeight*TankLinerThickness*TankCircumference*TankLinerDensity+8*TankLinerThickness*TankCrossSection*TankLinerDensity);
	//TankStructureMass = (TankHeight*ActualTankThickness*TankCircumference*TankDensity+8*ActualTankThickness*TankCrossSection*TankDensity+TankCircumference*2*TankRadius*MinimalMaterialThickness*TankDensity);//  cylindrical tanks with spherical tops and bottoms  with a cylinder connecting the two tanks     zzz add 1mm of aluminum to the insides of the tanks!
	TankStructureMass = (TankHeight*ActualTankThickness*TankCircumference*TankDensity+6*ActualTankThickness*TankCrossSection*TankDensity);//  cylindrical tanks with spherical tops and bottoms   and a common bulkhead

	//TankMass =  (TankStructureMass +RequiresLiner?TankLinerMass:0);
	TankMass =  TankStructureMass ;
        VehicleDryMass = (TankMass+EngineMass);
 /*  */
#define DesiredMaxAspectRatio 10. // basically disabled  we could also put a  angle of attack into the aerodynamica
//#define AspectWeightingFactor .0
#define AspectWeightingFactor .5
//aspectBogosityFactor = fmax(0.  ,1+TankHeight/(TankRadius *2)- DesiredMaxAspectRatio)*PayLoad*AspectWeightingFactor;
aspectBogosityFactor = fmax(0.  ,1+TankHeight/(TankRadius *2)- DesiredMaxAspectRatio)*PayLoad*AspectWeightingFactor;


	{
		ISParray [NSec]= 0;
		Xarray[NSec] = X;
		Yarray[NSec] = Y;
		XVelocityarray[NSec] = VelocityX;
		YVelocityarray[NSec] = VelocityY;
		Phiarray[NSec] = 0;
		Thrustarray[NSec] = 0;
	}
	while(VelocityY>=-1.0 && Y>-1.){// we let VelocityY go negative because we havent tested  centrifugal force yet NOte add that  
		if(MadeOrbit) {break;} // comment this out for full trajectory
		if(VelocityX <-FUZZ){break;}
		//if(Theta>(Pi/2)){break;}
		if(NSec >MaxSECS){break;}
		if(NSec>1200 && NSec >OrbitAltitude/500){break;}
		ISP=CalcISP(AirPressure,ChamberPress);
		if(ISP >MaxISP){
			MaxISPAltitude = Y;
			 MaxISP = ISP;
		}
		if(ISP <MinISP){
			 MinISP = ISP;
		}
		if(ISP <180.){ // this was -1 but if ISP falls below 100 give up
//if(NSec !=0){printf("Nsecs = %d ISP = %f  Y = %f  velicityY %f\n",NSec , ISP, Y,VelocityY);}
		//	return -100000.;
			break;
		}
		
		Acc = ThrustAndMassToAcc ;
		if(Acc >FUZZ){
			GBogosityFactor+=((( fmax(0,Acc -MaxDesiredAcc)+fmax(0,MinDesiredAcc -Acc))/GEarth)*PayLoad);  // integrate out of range acceleration and normalize wrt rocket size by factoring in payload
		}
		TotalDrag += Drag;
		sinOrbit = VelocityX/(EarthRadius +Y); /* small ang
les sin theta = theta  = VelocityX/radius  */
		yComponent =  - GEarth*(EarthRadius*EarthRadius)/((EarthRadius + Y)*(EarthRadius +Y)) + VelocityX*sinOrbit;
		VelocityY += sin(Theta)*Acc  + yComponent;
		xComponent = sinOrbit*VelocityY;

		VelocityX += cos(Theta)*Acc  - xComponent ;

/* here we need to transform vy & vy  */
		Velocity = sqrt((VelocityX-EarthRotationalVelocity)*(VelocityX-EarthRotationalVelocity)+VelocityY*VelocityY);
		Y += VelocityY;
		X += VelocityX;
		DeltaV += Acc;
		MFR=CalcMassFlow;
		VMass=VehicleMass;
		FPress=TankBottomPress;
		Thrust=CalcThrust;
		if(Thrust >MaxThrust){
			MaxThrust = Thrust;
		}
		THRUSTSUM += Thrust;
		if(NSec == 1){TakeOffThrust = Thrust;}

                       FractionOfDesiredHeightCalculated =(.5*VelocityY*VelocityY/GEarth  +Y ) /OrbitAltitude;





		if(Velocity<FlatVelocity){
                         /*Theta = (Pi/2)*(1-(Velocity)/(FlatVelocity));*/
                         Theta = (Pi/2)*(1-(((Velocity)/(FlatVelocity))*((Velocity)/(FlatVelocity))));
		}else{
		   Theta = 0;
		}


if(debugPrint&&Thrust >0){
printf("T %d X %G Y %G  A %G F %G Th %G VX %G VY %G Thru %G D %G\n",NSec,X,Y,Acc,FuelMass,Theta,VelocityX,VelocityY,Thrust,Drag);
}
		//Acc = ThrustAndMassToAcc ;
		FuelMass -= MFR;
		if(!Burnout &&(FuelMass <MFR || MFR == 0)) { OutOfFuel = 1;Burnout = NSec; BurnoutAltitude = Y;BurnoutVelocity = Velocity;}
		if(Velocity > MaxV) MaxV = Velocity;
		if(Y > MaxY) MaxY= Y;
		if(abs(Acc)>MaxAcc) MaxAcc = abs(Acc);
		if(abs(Acc)>.01 && abs(Acc)<MinAcc){
			 MinAcc = abs(Acc);
			MinAccAltitude = Y;
			MinAccVelocity = Velocity;	
		}
		if(FPress>MaxFPress) MaxFPress = FPress;
			if(Drag>MaxDrag) {
		  MaxDrag = Drag; 
		  MaxDragAlt = Y;
		MaxDragV = Velocity;
		}
		NSec +=1;
	{
		ISParray [NSec]= ISP;
		Xarray[NSec] = X;
		Yarray[NSec] = Y;
		XVelocityarray[NSec] = VelocityX;
		YVelocityarray[NSec] = VelocityY;
		Phiarray[NSec] =  RadiansToDegrees( /*Pi/2 -*/ Theta);
		Thrustarray[NSec] = Thrust;
	}
	}
    *cost_flag = TRUE;
//  if(Y<OrbitAltitude||VelocityX<OrbitalVelocityAtAltitude) *cost_flag = FALSE;
/*if(PascalsToPsi(MaxFPress)>=2000.0) { 

	*cost_flag = FALSE;
	Y = X  = VelocityX = VelocityY = 0;
	
}
*/
testNumber +=1;
if(debugPrint){
	printf("TankHeight %G TankCircumference %G MassFlowAllowed %G FlatVelocity %G  Merit %G\n",TankHeight,TankCircumference,MassFlowAllowed,FlatVelocity,Merit); 
}
//printf("testNumber = %d %G  \n",testNumber,Merit);
if(MaxMerit<Merit /*|| (Merit >0 && .9* MaxMerit<Merit)*/){ 
printedTestNumber +=1;
putchar(0x1b);
putchar(0x5b);
putchar(0x3b);
putchar(0x48);
putchar(0x1b);
putchar(0x5b);
putchar(0x32);
putchar(0x4a);
printf("\n\n");
	if(MaxMerit<Merit){
		MaxMerit=Merit;
/*printf("AN IMPROVEMENT\n");*/
	}else{
		printf("NOT BETTER BUT PERHAPS STILL INTERESTING %G  as good as best \n",100*Merit/MaxMerit);
	}

printf("Parameters from the simulated annealer for test# %d bestResult# %d:\n",testNumber,printedTestNumber);
printf("Expansion Ratio %G \n", ExpansionRatio);
printf("MaxDesiredG = %G \n",MaxDesiredG);
printf("TankHeight=  %Gm ",TankHeight);
printf("TankCircumference=  %Gm",TankCircumference);
printf(" TankRadius = %Gm",TankRadius);
printf("   aspect height/diameter %G\n",1+TankHeight/(TankRadius *2));
printf("Device Height %Gm %Gft+ engine and Cone  diameter = %Gm %Gft\n",  TankHeight+2*TankRadius,(TankHeight+2*TankRadius)*3.280840,TankRadius*2,TankRadius*2*3.280840);
printf("MassFloweAllowed %Gkg/s ", MassFlowAllowed);
printf("     calculated approx MaxThrust %G N  %G lb\n",MaxThrust,NewtonsToLbf(MaxThrust));

printf("   Flight path parameters ");
printf("      FlatVelocity=  %G\n",FlatVelocity);
printf("\n");
printf("VelocityX = %Gm/s Altitude= %Gm VelocityY %Gm/s Velocity %Gm/s\n",VelocityX,Y,VelocityY,Velocity);
printf("OrbitalVelocityAtAltitude %G ",OrbitalVelocityAtAltitude);
printf("OrbitAltitude %G m\n",OrbitAltitude);
printf("sumISP=  %G   MaxAcc=  %Gg MinAcc %Gg   MaxFPress=  %Gpsi\n",THRUSTSUM/(StartFuelMass*GEarth),MaxAcc/GEarth,MinAcc/GEarth,PascalsToPsi(MaxFPress));
//printf("MinAccAltitude %gm MinAccVelocity %Gm/sec \n",MinAccAltitude,MinAccVelocity);
printf("Burnout at %d seconds   sec velocity = %G altitude = %G modeling ends at %d sec\n",Burnout,BurnoutVelocity, BurnoutAltitude,NSec);

//printf("\nA fews constants\n");
printf("tankage pressure set to %G pascals %G psi ", TankPress, PascalsToPsi(TankPress));
printf("  PayLoad=  %Gkg\n",PayLoad);
printf("FuelDensity %G kg/(m^3) \n",FuelDensity);
printf("\ncost assumptions:  LaborFactor %G set to 1 for experimental FuelUnitCost $%G/m^3\n  EnginePrice $%G/kg  %sPrice $%G/kg\n",LaborFactor,FuelUnitCost,EngineCostPerWeight,MaterialName,MaterialPrice);
printf("\nTakeOffThrust  %Gn %G lb   Thrust/Weight   %G \n",TakeOffThrust,NewtonsToLbf(TakeOffThrust	),TakeOffThrust/(GEarth*(VehicleDryMass+PayLoad+StartFuelMass)));
printf("MaxThrust = %G n %G lb \n",MaxThrust, NewtonsToLbf(MaxThrust));
//printf("emperical thrust to weight ratio of engine = %G  we started with %G\n", (MaxThrust /GEarth	)/EngineMass,  ThrustWeightRatioOfEngine);

printf("TankCrossSection = %G\n",TankCrossSection);
printf("StartFuelMass   %Gkg  \n",StartFuelMass);
printf("TankMass %Gkg\n",TankMass);
printf("EngineMass %Gkg \n",EngineMass);
printf("\nTankHeight   %Gm TankRadius   %Gm TankVolume   %Gm^3\n",TankHeight,TankRadius,TankVolume);
if(EngineMass != NominalEngineMass){
	printf("    minimal engine used, calculated engine mass was:  NominalEngineMass %gkg\n",NominalEngineMass);
}
printf("ActualTankThickness       %G m  \n",ActualTankThickness);
//printf("AlternateTankThickness    %G m\n",AlternateTankThickness);
printf("\n");


printf("MaxISP = %G  ",MaxISP);
printf("MaxISPAltitude = %G  ",MaxISPAltitude);
printf("MinISP = %G  \n",MinISP);


//printf(" StartFuelMass  %G ResidualFuelMass %G MassRatio %G\n",StartFuelMass,FuelMass,(VehicleDryMass+PayLoad+StartFuelMass)/(VehicleDryMass+PayLoad));
//printf(" StartFuelMass  %G\n",StartFuelMass);


printf("DryMass  %Gkg DryMass+PayLoad %Gkg  PayLoad %Gkg StartFuelMass %Gkg   Vehicle_Loaded%G\n",VehicleDryMass, VehicleDryMass+PayLoad,PayLoad,StartFuelMass, VehicleDryMass+PayLoad+StartFuelMass);
printf("MassRatio  %G   DeltaV   %G  \n\n",(VehicleDryMass+PayLoad+StartFuelMass)/(VehicleDryMass+PayLoad), DeltaV);

	checkMR( THRUSTSUM/(StartFuelMass*GEarth),(VehicleDryMass+PayLoad+StartFuelMass)/(VehicleDryMass+PayLoad),DeltaV);
//printf("VehicleCost  $%G    Un//insuredLaunchCost  $%G\n",VehicleCost,UninsuredLaunchCost);
printf("\nVehicleCost  $%G    Tank $%G Engine $%G Fuel  $%G\n",VehicleCost, TankCost, EngineCost, FuelCost);
if (LaborFactor != 1.){printf("  MeritVehicleCost  with LaborFactor for the anealer %G\n",MeritVehicleCost);}
 printf(" NormalizedMaxDrag = %G MaxDrag=%G at Alt=%G m going %G M/s\n",MaxDrag/StartFuelMass,MaxDrag,MaxDragAlt,MaxDragV);
printf("GBogosityFactor = %G\n",GBogosityFactor);
printf("aspectBogosityFactor %G\n",aspectBogosityFactor);
printf("TotalDrag %G\n",TotalDrag);
printf("Profit= $%G\n",Profit);
printf("LaunchCost $%G\n\n\n",LaunchCost);

printf("MadeOrbit %s Merit %G MeritRatio %G %s  ",MadeOrbit?"Yup":"NO GODDAMNIT!",Merit,Merit/MaxMerit,Merit == MaxMerit?"BEST":"GOOD");
printf("aspect height/diameter %G\n\n",1+TankHeight/(TankRadius *2));
printf("Fudge = %G ",Fudge);
printf("engine thrust/weight ratio %G\n",ThrustWeightRatioOfEngine);
printf("MaxFPress %G ",PascalsToPsi(MaxFPress));
printf("TankPress %G ",PascalsToPsi(TankPress));
printf("MaxFPress/TankPress = %G \n",MaxFPress/TankPress);
printf("\n");
printf("a = %G b=%G c=%G d=%G \n",a,b,c,d);
printf("payload/(tank+engine) = %G inverse = %G",PayLoad/(TankMass+EngineMass),(TankMass+EngineMass)/PayLoad);
printf(" Latitude = %G ",Latitude);
printf("Chamber Pressure = %G \n",PascalsToPsi(ChamberPress));
printf("PseudoVolume = %G \n",TankPseudoVolume);
printf("$/pound =$ %G\n",LaunchCost/(PayLoad*2.2));
	for( i=0;i<MaxSECS&&i<=NSec;i++){
		if(Thrustarray[i] || !(i %1000)|| i+100>NSec){
//		printf("t = %i",i);
//		printf(" X = %G Y = %G XVel = %G YVel = %G ISP = %G Phi = %G Thrust = %G \n", Xarray[i], Yarray[i], XVelocityarray[i],YVelocityarray[i],ISParray[i], Phiarray[i], Thrustarray[i]);
	}
fflush(stdout);
}

fflush(stdout);
}
/* End meat of cost function. */
  funevals = funevals + 1;

#if TIME_CALC
  if ((PRINT_FREQUENCY > 0) && ((funevals % PRINT_FREQUENCY) == 0))
    {
      fprintf (ptr_out, "funevals = %ld  ", funevals);
      print_time ("", ptr_out);
    }
#endif
fflush(stdout);
  return (-Merit);
}

void checkMR(double isp, double massratio, double deltav){
printf(" As a cross check calculating delta v from ISP and Mass ratio\n");
printf("ISP %G  MassRatio %G  ln(MassRatio) * ISP *Gearth = %G\n",isp,massratio,log(massratio)*isp*GEarth );

}

static inline double calcExpandedPress(double, double, double );

 double internal_calcExpandedPress(double press, double expansionRatio, double ksp){
double guessIPress = 0;
static double oldguess = 0;
	double Uguess = press;
	double Lguess = press/expansionRatio;
        double delta;
	if(oldguess >Lguess && oldguess <Uguess){
		guessIPress = oldguess; // start with where we were last iteration
// and the pressure will rarely go up!
	}
	for(;;){

		delta = 1./expansionRatio - pow((ksp +1)/2,1/(ksp -1))*pow(guessIPress/press,1/ksp)*sqrt(((ksp+1)/(ksp-1)) * (1 - pow((guessIPress/press),(ksp-1)/ksp)) ); /* from sutton equ 3-25 p 55 */
		if(delta <0){
			Uguess = guessIPress;
		} else {
			Lguess = guessIPress;
		}
		if (Uguess -Lguess <.001){ 
	//printf("delta %f guess %f %f", delta, Lguess,Uguess);
	oldguess = Lguess;
			return Lguess;
		}
		guessIPress = (Lguess +Uguess)/2;
	}
}

double CF_ThrustCoefficient(double chamberpressure,double expandedpressure,double externalpressure,double ksp,double area_expansion_ratio){
  return  (sqrt((2.*ksp*ksp/(ksp-1.))*pow((2./(ksp+1.)),((ksp+1.)/(ksp-1.)))*(1.-pow((expandedpressure/chamberpressure),((ksp-1.)/ksp))))
        +((expandedpressure-externalpressure)/chamberpressure)*area_expansion_ratio); /* sutton equ 3-30 p 59 ?? */
}

double CalcISP (double pressure, double chamberPress){
static double oldpressure;
static double oldchamberPress;
static double oldresult;
double result;

  double expandedPressure;
  double thrustCofficient;
/* since I have coded this very inefficiently and it is inner loop and
it is called 4 times with the same arguments: I have built this simple cache.
If called with the identical arguments as last time, return the same results.
THis is almost safe, all the other parametere (yuck there are some)
are globals and should represent constatn aspects of a particular rocket flight 
*/
  if((pressure == oldpressure) && (chamberPress == oldchamberPress)){
	return oldresult;
   }else {
      oldpressure = pressure;
      oldchamberPress = chamberPress;
   }
 expandedPressure =  calcExpandedPress(chamberPress,ExpansionRatio,KSpecificHeatRatio) ;
    thrustCofficient =  CF_ThrustCoefficient(chamberPress,expandedPressure,pressure,KSpecificHeatRatio,ExpansionRatio);
    result =  CharisticISP * thrustCofficient;
//fprintf(stderr,"cf %G  isp %G exp %G cp %G externalpress %G\n",thrustCofficient,  CharisticISP * thrustCofficient,expandedPressure,chamberPress, pressure);
/*if(thrustCofficient <0 && NSec != 0){
printf("expanded press %f chamberPress %f CF = %f ISP = %f expansionRatio %f, externalPress %f NSec = %i \n", PascalsToPsi(expandedPressure),PascalsToPsi(chamberPress),thrustCofficient, result, ExpansionRatio, pressure, NSec);
}*/
  oldresult = result*Fudge;
  return result *Fudge;
}

/*
double CalcISP3 (double pressurePascal, double chamberPress){ // hard coded for 3000 psi combustion chamber 
double r = 20;
double p = 2500;
double result;
double pressure;
	pressure = PascalsToPsi (pressurePascal);
	if (ExpansionRatio>12.159){
	p = 15; r = 3037.7;
        }
	if (ExpansionRatio>16.496){
        p= 10  ;r= 3103.45 ;
        }
	if (ExpansionRatio>27.851){
        p=5 ;r  = 3201.71 ;
        }
	if (ExpansionRatio>32.977){
        p=4   ;r= 3229.9 ;
        }
	if (ExpansionRatio> 41.005){
        p=3   ;r= 3264.23 ;
        }
	if (ExpansionRatio>55.747){
        p=2  ;r = 3308.61;
        }
	if (ExpansionRatio> 5.7006*12.159){
        p=1.5 ;r= 3337.54;
        }
	if (ExpansionRatio>12.159*7.7453){
        p=1  ;r = 3375.04;
        }
	if (ExpansionRatio>12.159* 13.059){
        p=.5 ;r= 3431.13;
        }
	if (ExpansionRatio>12.159* 25.955){
        p=.2  ;r= 3491.89 ;
        }
	if (ExpansionRatio>12.159* 43.517){
        p=.1  ;r= 3529.38;
        }
	if (ExpansionRatio>12.159* 43.517* 1.6724){
        p=.05 ;r= 3560.77;
        }
	if (ExpansionRatio>12.159* 43.517*3.2875){
        p=.02 ;r= 3594.45;
        }
	if (ExpansionRatio>12.159* 43.517* 5.4638){
        p=.01 ;r= 3615.03 ;
        }
	if (ExpansionRatio>12.159* 43.517* 5.4638* 5.2902){
        p=.001 ;r= 3665.17 ;
        }
	if (ExpansionRatio>12.159* 43.517* 5.4638* 407270.){
        p=.00000000008 ;r= 3708.78 ;
        }
//printf("pressure = %G p = %G foo = %G\n",pressure,p,((p-pressure)) * ExpansionRatio);// is the following close to the right units?
#define SpecificArea .000013 // area of throat for 1kg propelant!
    result = (r * Fudge * (3135.35/3037.7) + PsiToPascals(p-pressure)*ExpansionRatio*SpecificArea)/GEarth;  
    //return (r * Fudge * (3135.35/3037.7) )/GEarth;  
//printf("pressure =%f", PascalsToPsi(chamberPress));
//printf("approx G throttle %f ", chamberPress/ChamberPressConstant);
	if(result <0.){
		return 0.;
	}else {
		return result;
	}
}*/
// fix to smooth it out 
// add water enhancement
// trajectory adjustment stuff
static inline double Throttle(double y, double a1 ,double b1, double c1,double d1)
{
// throttle depend on altitude !!
return 1.; // commented out 
if(y<a1){
	return 1;
}
if( y < a1 + b1){
	return fmin(1. , c1);
}
if( y < a1 + b1+b1){
	return fmin(1.,d1);
}
  return 1.;
//return  fmax(0., fmin (1., a1 * y*y*y + b1*y*y + c1*y + d1));

}
/* ***********************************************************
// The caclulation above, of expanded pressure is quick and dirty, and probably
// too precise, but I'll not bother to spped it up  algorithmically
// instead I'm building a table of results for a constatn Kspecificheat
// on the fly, then we just pick the entry in the table without even bothering
// to interpolate.  the table is pretty flat except near zero pressure and 
// expansion ratio.  I expect that those results won't get to orbit, but they 
// might contribute to hill climbing appropriately.
// Nope I changed my mind (rusults!!)  I'll interpolate.  the asa
// seems to find the nonlinearity !!!
************************************************************** */
#define NExpansionBuckets 250
#define NCombustionPressureBuckets 100
#define MaxExpansionRatio 1000
#define PressureIncrement (ChamberPressConstant/NCombustionPressureBuckets)
#define ExpIncrement (MaxExpansionRatio/NExpansionBuckets)
double ExpPressArray[NExpansionBuckets][NCombustionPressureBuckets];

void setupExpArray(){
int i,j;
     for(j = 0;j<NCombustionPressureBuckets;j++){
	ExpPressArray[0][j] = 0.;
     }
     for(i = 1;i<NExpansionBuckets;i++){
//printf("expansion = %d %d \n",i*ExpIncrement,i);
	ExpPressArray[0][j] = 0.;
	for(j = 1;j<NCombustionPressureBuckets;j++){
	    ExpPressArray[i][j] =  internal_calcExpandedPress(j*PressureIncrement, i*ExpIncrement,KSpecificHeatRatio);
//printf(" %d %f ", j, PascalsToPsi(ExpPressArray[i][j]));
	}
//printf("\n");
    }
}

static inline  double calcExpandedPress(double press, double expansionRatio, double ksp){
int expIndex;
int pIndex;
double ret;
double r;
        expIndex = imin((int)(expansionRatio/ExpIncrement),NExpansionBuckets-1);
	pIndex   = (int)(press/PressureIncrement);
	r = (press/PressureIncrement) - pIndex;
	ret =  ExpPressArray[expIndex][pIndex]
	 + (ExpPressArray[expIndex+1][pIndex+1]
	- ExpPressArray[expIndex][pIndex]) * r;
	//- ExpPressArray[expIndex][pIndex]) * ((((int)expansionRatio)%((int)ExpIncrement))/ExpIncrement);

//fprintf(stderr,"press = %G exratio = %G cf = %G expIndex %d pIndex %d  exparray = %G r = %G\n",press , expansionRatio,ret,expIndex, pIndex,ExpPressArray[expIndex][pIndex],r);
	return ret;
}