// major rev, sept 4 2001 changed  hillclimbing variables to orthagonalize better
//mostly changet it to work directly from volume and aspect ratio, rather
//than height and width, thus the anealer doesent have to autocorelate these
//which it doesn't really, factor of 10 improvement in convergence, this
// allows other nonconvergent off the wall things to be attempted, like
// like the thrust and coast mechanism that was in a previous version.

/* TODO:
//  simplyfy pressure calculation
//  rationalize the things that the asa varries eg) volume, thrust, aspect-ratio,tank/engine-dryweight-ratio,takeoff-g
// output a sepreate file of vrml models!  this is trivial, so do it!
// clean up or remove the pressure with acc stuff, pressure may increase slightly with increase in isp over takeoff isp that is all so simplify it
// fix pressure stuff to require 
*/
#include "asa_user.h"
#include <math.h>
#include <stdio.h>
double payload;
#define PayLoad payload //shared with user.c !!
#define AluminumDensity (2702.0) //C//CONSTANTonstant
#define PET 0
#define PETDensity       (1370.)////CONSTANTConstant
#define PETStrength       (3.44379e8)////CONSTANTConstant
#define PETPrice	   (5.) // Constant guess
#define ALUMINUM 2
#define ALUMINUMDensity  (2702.0) //C//CONSTANTonstant
#define ALUMINUMStrength  (6.e8)  //random guess enhanced for  cryotemp
#define ALUMINUMPrice	 (5.)    // another guess	
#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 MATERIAL PET
//#define MATERIAL ALUMINUM
//#define PF //PRESSUREFED
#ifdef PF //PRESSUREFED
#define CHAMBERPRESSURE (TankPressureexternal *.8)
#define PRESSUREIZATIONREQUIRED 13
#else
#define CHAMBERPRESSURE 10000
#define PRESSUREIZATIONREQUIRED 20
#endif
#define ChamberPressConstant ((PsiToPascals(CHAMBERPRESSURE)))//Parameter
#define  CharisticThrustWeightRatioEngine (220.) // 100 high normal range we might be 300 or 1000//Parameter
#define Fudge .9//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//need more data  from proprep especially for pressure fed  stuff

#define KSpecificHeatRatio (1.24)  //Parameter sutton table 6-4 p 194
// KSpecificHeatRatio is the specific reat ratio as defined by KSpecificHeadRatio
 #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);
static double oldpressure = -1;
static double oldchamberPress = -1;
int ExpArrayInitialized = 0;
//#define MaxSECS 140002
#define MaxSECS 1000000
#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 FUZZ .0001//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 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

//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 )
#define CalcThrust (CalcMassFlow*CalcISP(AirPressure(Y),ChamberPress)*GEarth)
//#define CalcMassFlow fmin(CharisticThrustWeightRatioEngine*EngineMass/CharisticISP,(OutOfFuel?0.:(TankCrossSection/ExpansionRatio)*(ChamberPress/CalcISP(AirPressure(Y),ChamberPress))))// stop burning when out of fuel(AirPressure(Y))
//#define CalcMassFlow fmin(CharisticThrustWeightRatioEngine*EngineMass/CharisticISP,(OutOfFuel?0.:(TankCrossSection/ExpansionRatio)*(ChamberPress/CalcISP(AirPressure(Y),ChamberPress))))// stop burning when out of fuel(AirPressure(Y))
//#define CalcMassFlow ((fmin(CharisticThrustWeightRatioEngine*EngineMass/CharisticISP,(OutOfFuel?0.:1000000. /* sure to be big*/))),printf("MassFlow = %g OutOfFuel %i\n",CharisticThrustWeightRatioEngine*EngineMass/CharisticISP,OutOfFuel))// stop burning when out of fuel(AirPressure(Y))
//#define CalcMassFlow (printf("MassFlow = %g OutOfFuel %i\n",CharisticThrustWeightRatioEngine*EngineMass/CharisticISP,OutOfFuel),(fmin(CharisticThrustWeightRatioEngine*EngineMass/CharisticISP,(OutOfFuel?0.:1000000. /* sure to be big*/))))
#define CalcMassFlow (fmin(CharisticThrustWeightRatioEngine*EngineMass/CharisticISP,(OutOfFuel?0.:1000000. /* sure to be big*/)))
//#define CalcMassFlow (fmin(CharisticThrustWeightRatioEngine*EngineMass/CharisticISP,(OutOfFuel?0.:1000000. /* sure to be big*/)))// stop burning when out of fuel(AirPressure(Y))

#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* (1-PressureTankFraction))
//#define FuelCost (FuelUnitCost*TankVolume)
#define FuelHeight  (FuelVolume/TankCrossSection)
//#define TankBottomPress (TankPress+Acc*(FuelHeight*FuelDensity))//do we have teh units right here? zzz

#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 double
	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;
}
// IMiles and IOrbit are to produce an integer for the #if
#define MilesToM(x) (x*1609.344)
#define IMilesToM(x) (x*1609)//CONSTANT integer expression for defines!
//#define MilesOrbitalAltitude 25000
#define MilesOrbitalAltitude 100
#define OrbitAltitude             MilesToM(MilesOrbitalAltitude)//Parameter // lunar!
#define IOrbitAltitude             IMilesToM(MilesOrbitalAltitude)//Parameter // lunar!
//#define OrbitAltitude  200000.  
//#define OrbitAltitude  100000.  // X prize altitude
//#define OrbitAltitude  200000.  
//#define OrbitAltitude  400000.  // sounding rocket CAT prize altitude
#define GeosyncAltitude             35786000  //geosync//Parameter
#if (IOrbitAltitude >(GeosyncAltitude +10))
#define AltitudeOnly  TRUE // only try to match altitude, current design of trajedctory is crude, this helps.
#else // pick one of the following according to taste
//#define AltitudeOnly  TRUE
#define AltitudeOnly  FALSE
#endif
#define MadeOrbit  ((AltitudeOnly?1:(VelocityX>OrbitalVelocityAtAltitude))&&(Y>OrbitAltitude))
#define MadeAltitude  ((Y>OrbitAltitude))
#define Price (MadeOrbit?PayLoad*(2000.0*2.2):0)////Parameter
#define Profit (Price-LaunchCost)
/* Merit is used by the asa to pick better vehicles to test
it is altitude till it makes orbit or altitude
then is cost perhaps midified by residual y velocity
the y velocity factor is probably only  needed because of the semi bogus attidude policy 
that the code currently uses tradeoff seems to be at 100 grams  the real fix is to make a better
attidude plan TODO zzz
*/
//#define Yvelocitybogusfactor (PayLoad <.11 ?0:VelocityY/100) 
#define Yvelocitybogusfactor 0
#if AltitudeOnly 
#define Merit ((MadeOrbit?(Profit/UninsuredLaunchCost-Yvelocitybogusfactor)-aspectBogosityFactor:(fmin(OrbitAltitude,fmax(1,Y))-10e5)))   //
#else
//#define Merit (MadeOrbit?(Profit/UninsuredLaunchCost- Yvelocitybogusfactor+1000):(((fmax(0,fmin(1,(Y/OrbitAltitude)))*1000 +fmax(0,fmin(1,VelocityX/OrbitalVelocityAtAltitude))*1000-11e6))))//
//#define Merit (MadeOrbit?(Profit/UninsuredLaunchCost- Yvelocitybogusfactor+1000):((MadeAltitude)?((fmin(OrbitAltitude,fmax(1,Y))/20 +fmin(OrbitalVelocityAtAltitude,fmax(1,VelocityX))-11e6)):(fmax((VelocityX),1)-7e6)))//
//#define Merit (MadeOrbit?(Profit/UninsuredLaunchCost- Yvelocitybogusfactor+1000):((!MadeAltitude)?((fmin(OrbitAltitude,fmax(1,Y))/.2000 +fmin(OrbitalVelocityAtAltitude,fmax(1,VelocityX))-11e6)):(fmax((VelocityX),1)-7e6)))//
//#define Merit (MadeOrbit?(Profit/UninsuredLaunchCost- Yvelocitybogusfactor):(fmin(OrbitAltitude,fmax(1,Y)) +fmin(OrbitalVelocityAtAltitude,fmax(1,VelocityX))-7e5)) //
#define Merit (MadeOrbit?(Profit/UninsuredLaunchCost- Yvelocitybogusfactor):(fmin(OrbitAltitude,fmax(1,Y))/10 +fmin(OrbitalVelocityAtAltitude,fmax(1,VelocityX))-7e5)) //

#endif

#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 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

/*
double Compute(double h) {
// this foolish thing calculates in English Units  but i'll leave it that way for now

   double v = 0;
//var rl = document.forms[1].rlength.value;
double TEMPSL = 518.67;
//double RHOSL = 0.00237689;
//double PRESSSL = 2116.22;
double saTheta = 1.0;
double saSigma = 1.0;
double saDelta = 1.0;
if ( h<232940 ){
saTheta = 1.434843 - h/337634;
saSigma = pow( 0.79899-h/606330, 11.20114 );
saDelta = pow( 0.838263-h/577922, 12.20114 );
}
if ( h<167323 ){
saTheta = 0.939268;
saSigma = 0.00116533 * exp( (h-154200)/-25992 );
saDelta = 0.00109456 * exp( (h-154200)/-25992 );
}
if ( h<154199 ){
saTheta = 0.482561 + h/337634;
saSigma = pow( 0.857003+h/190115, -13.20114 );
saDelta = pow( 0.898309+h/181373, -12.20114 );
}
if ( h<104987 ){
saTheta = 0.682457 + h/945374;
saSigma = pow( 0.978261+h/659515, -35.16319 );
saDelta = pow( 0.988626+h/652600, -34.16319 );
}
if ( h<65617 ){
saTheta = 0.751865;
saSigma = 0.297076 * exp( (36089-h)/20806 );
saDelta = 0.223361 * exp( (36089-h)/20806 );
}
if ( h<36089 ){
saTheta = 1.0 - h/145442;
saSigma = pow( 1.0-h/145442, 4.255876 );
saDelta = pow( 1.0-h/145442, 5.255876 );
}
{
  double tempVal = TEMPSL * saTheta;
//double rhoVal = RHOSL * saSigma;
//double pVal = PRESSSL * saDelta;
//double viscVal = 0.0226968*pow( tempVal, 1.5 ) / ((tempVal)+198.72) /1000000.0;
double soundVal = sqrt( 1.4*1716.56*(tempVal) );
double machVal = v/soundVal;
//double qVal = 0.7*pVal*machVal*machVal;
//double rlength = 1.0; // zzz hack may need dimension fixes
//double reynolds = v*rlength*rhoVal/viscVal;

//double cfturb = 0.455/pow((log(reynolds)/log(10)),2.58);

//double cflam = 1.328/sqrt(reynolds);

double cpmin;
//double cpstar =  (pow((1/1.2 +machVal*machVal/6.0),3.5)-1)/(0.7*machVal*machVal);

cpmin =  -1.0/(0.7*machVal*machVal);

}
 return 0;
}
*/

  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.
nt)//
*****************************************************************/
static /*inline*/ double air_(double y){  // this takes 30% of the cpu!
  int i;
  double a = 0.;
  
  int km;
	//i = y/1000.;
   if(y>30000){// if    we are way up approximate  more
	if(y>200000.){
		return 0.;
	}
	i  = y/1000;
	return airdens[i]/2.2;
   }
	if (y<0){
		return airdens[0]/2.2;
		//i = 0;
	}	
	//i = ((int)y)/1000;
	i = y/1000;
//	km = y/1000.; // hack it to do more in int!
//return airdens[i];
	//km = y;
	//km = km/1000;
	km = i; // because we had it lying around
	a =  (y - km * 1000)/1000.;
	  return (airdens[i] - a*(airdens[i]  - airdens[i+1]))/2.2;
}

static double oldapy= -1;
static double oldapp= 0;
//airdens[0] = 1.225;
static /*inline*/ double AirPressure(double y) {
	if(y == oldapy){
		return oldapp;
	}
       //oldapp = PsiToPascals((air_(y)/airdens[0])* 15.);
       oldapp = PsiToPascals((air_(y)*2.2/1.225)* 15.);
	return oldapp;
}

extern int _flsbuf (unsigned int, FILE *);

  int NSec = 0;


double RadiusFrom(double volume,double aspectratio){
  double r;
	//calculate height from aspect ratio and volume
	// gemetry -s 2 spheres and a clyinder
	// vol = 2*4/3 *Pi *r^2 + Pi*r^2 * h
	// r = h/aspectratio
	// or h = r*aspectratio (simpler to write the solution this way)
	// vol = Pi*r^2(8/3r + r*aspectratio)
	r =  pow (volume/(Pi*aspectratio+6/3),.333333);
	return r;
}
double HeightFrom(double volume,double aspectratio){
  double r;
	//calculate height from aspect ratio and volume
	// gemetry -s 2 spheres and a clyinder
	// vol = 2*4/3 *Pi *r^2 + Pi*r^2 * h
	// r = h/aspectratio
	// or h = r*aspectratio (simpler to write the solution this way)
	// vol = Pi*r^2(8/3r + r*aspectratio)
	r =  pow (volume/(Pi*aspectratio+6/3),.333333);
	return 2*r + r*aspectratio;
}
static double ExpansionRatio;
static double TankPressureexternal;
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;
//proposed new set of x[]
double TakeoffG = x[0]; // perhaps not used
double TankVolume = x[1]*(1+( PayLoad/1000));//is as good as massratio ,but excludes otherdryweights
double etweightratio = x[2];
double AspectRatio = x[3];
  //ExpansionRatio = x[4];  
	double FlatVelocity = x[5];
#ifdef PF //PRESSUREFED //defined above
	double TankPressE = x[6];
#else
	double TankPressE = 0;
#endif
	TankPressureexternal = TankPressE;
	double TankPress = PsiToPascals(TankPressE+PRESSUREIZATIONREQUIRED);// for zzz now, there is nothing in the model to make the tank pressure stay high, note that we will have to multiplyin  an acceleration factor of MAXISP/CharisticISP but this may be quite close to 1
  int MaterialType = MATERIAL;//0 PET 1 Dyneema + liner ALUMINUM 2
double DeviceHeight = HeightFrom(TankVolume,AspectRatio);
double TankRadius = RadiusFrom(TankVolume,AspectRatio);
double  TankCrossSection = Pi*(TankRadius*TankRadius);
double CylinderHeight = 2*TankRadius*AspectRatio;
  double TankCircumference = 2* Pi * TankRadius;
//double TankCrossSection =  ((TankRadius*TankRadius)*Pi);
/* 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
// Initial conditions
//  double CylinderHeight = TankCircumference if CylinderHeight<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 = AltitudeOnly?0.:  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];
	double Drag = 0;

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 (AltitudeOnly){
FlatVelocity = 0; // just to show it isn't used in the calculations
}
	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 ALUMINUM:
		  MaterialPrice =ALUMINUMPrice;	
                  MaterialStrength = ALUMINUMStrength;
                  MaterialDensity = ALUMINUMDensity;
                  MaterialName = "ALUMINUM";
		  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;
	TankLinerMass =(CylinderHeight*TankLinerThickness*TankCircumference*TankLinerDensity+6*TankLinerThickness*TankCrossSection*TankLinerDensity);
	TankStructureMass = (CylinderHeight*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 ;
#ifdef PF
        //EngineMass = (TankMass *etweightratio); 
        EngineMass = (TankMass + StartFuelMass)*etweightratio; 
        //EngineMass = (StartFuelMass*etweightratio); 
        //EngineMass = etweightratio; 
        //EngineMass = 30; 
//printf("PF enginemass = %G etweightratio = %g tankmass = %g\n",EngineMass, etweightratio, TankMass);
#else
        EngineMass = TankMass*etweightratio;
#endif
        VehicleDryMass = (TankMass+EngineMass);
 /*  */
//#define DesiredMaxAspectRatio 10. // basically disabled  we could also put a  angle of attack into the aerodynamica
#define DesiredMaxAspectRatio 4. //
//#define AspectWeightingFactor .0
#define AspectWeightingFactor .5
aspectBogosityFactor = 0;


	{
		ISParray [NSec]= 0;
		Xarray[NSec] = X;
		Yarray[NSec] = Y;
		XVelocityarray[NSec] = VelocityX;
		YVelocityarray[NSec] = VelocityY;
		Phiarray[NSec] = 0;
		Thrustarray[NSec] = 0;
	}
	 //if (Profit/UninsuredLaunchCost <MaxMerit/1.1  && MaxMerit > 0) return -Profit/UninsuredLaunchCost;  /* don't bother*/
	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(Burnout && NSec <2){break;  };  //must be a bug!ZZZ
		if((!AltitudeOnly)&&VelocityX <-FUZZ){break;}
		//if(Theta>(Pi/2)){break;}
		if(NSec >MaxSECS){break;}
		if(NSec>1200 && NSec >OrbitAltitude/500){break;}
		ISP=CalcISP(AirPressure(Y),ChamberPress);
		if(ISP >MaxISP){
			MaxISPAltitude = Y;
			 MaxISP = ISP;
		}
		if(ISP <MinISP){
			 MinISP = ISP;
		}
		if(ISP <100.){ // 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);}
			break;
		}
		
		Drag = (CoefficientOfDrag/2.0)*air_(Y)*(Velocity)*(Velocity)*TankCrossSection ;// what units is Drag ? zzz  newtons??
		TotalDrag += Drag;
		Acc = ThrustAndMassToAcc ; // uses Drag 
		if(Acc >FUZZ){
			GBogosityFactor= 0;
		}
		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 AltitudeOnly
   Theta = Pi/2.0;
#else

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

		}
#endif

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;
	}
	}
if(NSec <2){
 Y =0; // break the trick of burning more fuel in 1 sec than we  have
Velocity = 0;
//aw fuck it just return a bad evaluation
//return -12e10;
}
    *cost_flag = TRUE;
//  if(Y<OrbitAltitude||VelocityX<OrbitalVelocityAtAltitude) *cost_flag = FALSE;
testNumber +=1;
if(debugPrint){
	printf("CylinderHeight %G TankCircumference %G  FlatVelocity %G  Merit %G\n",CylinderHeight,TankCircumference,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("VEHICLE DIMENSIONS\n\n");


printf("Expansion Ratio %G \n", ExpansionRatio);
//nprintf("MaxDesiredG = %G \n",MaxDesiredG);
printf("CylinderHeight=  %Gm \n",CylinderHeight);
printf("AspectRatio = %G \n",AspectRatio);
printf("TankCircumference=  %Gm",TankCircumference);
printf(" TankRadius = %Gm",TankRadius);
printf("\naspect height/diameter %G\n",1+CylinderHeight/(TankRadius *2));
printf("Device Height %Gm %Gft+ engine and Cone  diameter = %Gm %Gft\n",  CylinderHeight+2*TankRadius,(CylinderHeight+2*TankRadius)*3.280840,TankRadius*2,TankRadius*2*3.280840);
//printf("\ncalculated approx MaxThrust %G N  %G lb\n",MaxThrust,NewtonsToLbf(MaxThrust));
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  ",StartFuelMass);
printf("TankMass %Gkg ",TankMass);
printf("EngineMass %Gkg ",EngineMass);
printf("\nCylinderHeight   %Gm TankRadius   %Gm TankVolume   %Gm^3\n",CylinderHeight,TankRadius,TankVolume);
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  \nVehicle_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("tankage pressure set to %G pascals %G psi ", TankPress, PascalsToPsi(TankPress));
printf("tankpressure input = %G psi \n",TankPressureexternal);
printf("  PayLoad=  %Gkg\n",PayLoad);
printf("FuelDensity %G kg/(m^3) \n",FuelDensity);
//printf("aspect height/diameter %G\n\n",1+CylinderHeight/(TankRadius *2));
printf("Fudge = %G ",Fudge);
printf("charistic_engine thrust/weight ratio %G\n",CharisticThrustWeightRatioEngine);
printf("Max engine thrust/weight ratio %G\n",Fudge*CharisticThrustWeightRatioEngine*MaxISP/CharisticISP);
//printf("MaxFPress %G ",PascalsToPsi(MaxFPress));
printf("TankPress %G ",PascalsToPsi(TankPress));
//printf("MaxFPress/TankPress = %G \n",MaxFPress/TankPress);
printf("\n");
printf("Chamber Pressure = %G psi \n",PascalsToPsi(ChamberPress));
printf("Volume = %G \n",TankVolume);




printf("\n\nFlight parameters \n\n");
printf("      FlatVelocity=  %G\n",FlatVelocity);
printf("\n");
printf(" Latitude = %G ",Latitude);
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   \n",THRUSTSUM/(StartFuelMass*GEarth),MaxAcc/GEarth,MinAcc/GEarth);
if(!Burnout){printf("Fuel residual mass = %g\n", FuelMass);}
//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(" NormalizedMaxDrag = %G MaxDrag=%G at Alt=%G m going %G M/s\n",MaxDrag/StartFuelMass,MaxDrag,MaxDragAlt,MaxDragV);
//printf("GBogosityFactor = %G ",GBogosityFactor);
//printf("aspectBogosityFactor %G ",aspectBogosityFactor);
printf("TotalDrag %G\n ",TotalDrag);

//printf("\nA fews constants\n");
printf("\ncost assumptions:  LaborFactor %G FuelUnitCost $%G/m^3\n  EnginePrice $%G/kg  %sPrice $%G/kg\n",LaborFactor,FuelUnitCost,EngineCostPerWeight,MaterialName,MaterialPrice);


//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("Profit= $%G\n",Profit);
printf("LaunchCost $%G\n\n",LaunchCost);

printf("MadeOrbit %s Merit %G MeritRatio %G %s  ",MadeOrbit?"Yup":"NO GODDAMNIT!",Merit,Merit/MaxMerit,Merit == MaxMerit?"BEST":"GOOD");
printf(" MassRatio %G\n",(VehicleDryMass+PayLoad+StartFuelMass)/(VehicleDryMass+PayLoad));
printf(" PayloadRatio %G\n",(VehicleDryMass+PayLoad+StartFuelMass)/(PayLoad));
printf("etratio = %G\n",etweightratio);
//printf(" TakeoffG = %G\n",TakeoffG);
printf("TankVolumeabstract = %G\n",x[1]);
printf("ExpansionRatio = %G\n", ExpansionRatio);

//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("$/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);
oldpressure = -1;
oldchamberPress = -1;
//printf("Merit = %g \n",Merit);
  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;
	int i;
//	if(oldguess >Lguess && oldguess <Uguess){
//		guessIPress = oldguess; // start with where we were last iteration
// and the pressure will rarely go up!
//	}
	for(i=1;i++;){

		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){ // or .1 use this if we get the table driven form   working again
		//if (Uguess -Lguess <1.){ // actually probably good enough
	//printf("internal_calcExpandedPress delta %f guess %f %f %f i = %d\n", delta, Lguess,Uguess,Uguess -Lguess,i);
	oldguess = Lguess;
			return Lguess;
		}
		guessIPress = (Lguess +Uguess)/2;
	}
}

double CF_ThrustCoefficient(double chamberpressure,double expandedpressure,double externalpressure,double ksp,double area_expansion_ratio){
double thrustCofficient;
 thrustCofficient = (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 ?? */
if(thrustCofficient >5){
	printf(" thrustCofficient %f \n",thrustCofficient);
	printf(" chamberpressure %f expandedpressure %f externalpressure %f ksp %f area_expansion_ratio %f \n",  chamberpressure, expandedpressure, externalpressure, ksp ,area_expansion_ratio);
	printf("this is way too good, I consider it an error!;");
	//abort();
	return 0;
}
//	printf(" thrustCofficient %f \n",thrustCofficient);
//	printf(" chamberpressure %f expandedpressure %f externalpressure %f ksp %f area_expansion_ratio %f \n",  PascalsToPsi(chamberpressure), PascalsToPsi(expandedpressure), PascalsToPsi(externalpressure), ksp ,area_expansion_ratio);//zzz
return  thrustCofficient;
}

double CalcISP (double pressure, double chamberPress){
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;
   }
#ifdef PFx
 expandedPressure =  internal_calcExpandedPress(chamberPress,ExpansionRatio,(double)KSpecificHeatRatio) ;
#else
 expandedPressure =  calcExpandedPress(chamberPress,ExpansionRatio,(double)KSpecificHeatRatio) ;
#endif
    thrustCofficient =  CF_ThrustCoefficient(chamberPress,expandedPressure,pressure,KSpecificHeatRatio,ExpansionRatio);
    result =  CharisticISP * thrustCofficient;
if(thrustCofficient <1.3){
//fprintf(stderr,"cf %G  isp %G exp %G cp %G externalpress %G\n",thrustCofficient,  CharisticISP * thrustCofficient,expandedPressure,chamberPress, pressure);
}
if(thrustCofficient >5){
	printf("this is way too good, I consider it an error!;");
	abort();
}
/*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+1][NCombustionPressureBuckets+1];

void setupExpArray(){
int i,j;
//printf("setupExpArray \n");
     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;
 ksp = KSpecificHeatRatio;
        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,"calcExpandedPressZZZ press = %G exratio = %G cf = %G expIndex %d pIndex %d  exparray = %G r = %G result = %G \n",press , expansionRatio,ret,expIndex, pIndex,ExpPressArray[expIndex][pIndex],r,ret);
	return ret;
}