//// 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:
//add ullage and expecially pressurization gas calculations
// add staging
//  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!
*/
/* for 2 stage 
fix all wht weight and thrust calculations and acc to match new  conditions and add printouts.
add ullage and presurization weights!
 */
#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 ALUMINUMStrength  (3.e8)	//conservative 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 CHAMBERPRESSURE 3000
#define PRESSUREIZATIONREQUIRED 20
#endif
#define ChamberPressConstant ((PsiToPascals(CHAMBERPRESSURE)))	//Parameter

#ifdef PF
#define  CharisticThrustWeightRatioEngine (100.)	// 100 high normal range for pressure fed
#else
#define  CharisticThrustWeightRatioEngine (220.)	// 100 high normal range we might be 300 or 1000//Parameter
#endif



#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 ChamberPress (ChamberPressConstant )
#define CalcThrust (CalcMassFlow*CalcISP(AirPressure(Y),ChamberPress)*GEarth)
#define CalcMassFlow (fmin(CharisticThrustWeightRatioEngine*EngineMass/CharisticISP,(OutOfFuel?0.:1000000. /* sure to be big*/)))

#define ThrustAndMassToAcc ((CalcThrust-Drag)/CurrentVehicleMass)
#define TotalFuelVolume (TotalFuelMass/FuelDensity)
#define FuelUnitCost 500.0	// about $500/cubic meter of LOX (or fuel)//Parameter/CONSTANTramet///er//////
#define TotalTankVolume  (stage1TankVolume +stage2TankVolume)
#define   TotalFuelMass  (stage1FuelMass + stage2FuelMass)
#define TotalFuelCost (FuelUnitCost*TotalFuelVolume)
#define FuelHeight  (FuelVolume/TankCrossSection)
//#define TankBottomPress (TankPress+Acc*(FuelHeight*FuelDensity))//do we have teh units right here? zzz1GGgg
#define TotalTankMass (stage1TankMass +stage2TankMass)
//#define TotalVehicleMass  (TotalTankMass+TotalFuelMass+PayLoad)
#define TotalVehicleMass  (TotalTankMass+TotalFuelMass+PayLoad+ stagingPenalty)


#define stage2VehicleMass (stage2TankMass + stage2FuelMass + PayLoad)
#define MeritVehicleCost (TotalFuelCost+LaborFactor*(TotalTankCost+TotalEngineCost))
#define VehicleCost (TotalFuelCost+(TotalTankCost+TotalEngineCost))
#define stage1TankCost (stage1TankMass*MaterialPrice)
#define TotalTankCost (stage1TankCost +stage2TankCost)
#define EngineCostPerWeight 50.	//Parameter

#define stage1EngineCost  (fmax(EngineCostPerWeight*0.,(stage1EngineMass *EngineCostPerWeight)))

#define TotalEngineCost (stage1EngineCost + stage2EngineCost)
#define TotalEngineMass (stage1EngineMass + stage2EngineMass)
#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 t/ 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


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 + 8 / 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/3)r + r*aspectratio)
  r = pow (volume / (Pi * aspectratio + 8 / 3), .333333);
  return 2 * r + r * aspectratio;
}
static double ExpansionRatio;
static double TankPressureexternal;

double DeviceHeight;
double CurrentVehicleMass;
double TankRadius;
double TankCrossSection;
double CylinderHeight;
double TankCircumference;
double stage1TankMass;
double stage2TankMass;
double stageTankVolume;
double stage2FuelMass;
double stageFuelMass;
double stage1FuelMass;
double stage1EngineMass;
double stage2EngineMass;
double stage2TankVolume;
double stage1CylinderHeight;
double stage1TankRadius;
double stage2CylinderHeight;
double stage2TankRadius;

double stage1TankThickness;
double stage2TankThickness;
double stage1ExpansionRatio;
double stage2ExpansionRatio;

double stagingPenalty = 0;
double printableStagingPenalty = 0;
#define NStages 2

  double MaterialDensity = 0;

  int stage;
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 stage1TankVolume = x[1] * (1 + (PayLoad / 1000));	//is as good as massratio ,but excludes otherdryweights
  double stage1etweightratio = x[2];
  double stage1AspectRatio = x[3];
  double stage2etweightratio = x[8];
  double stage2AspectRatio = x[9];
  //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


/* 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;
  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.91e17;
  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;
  double MaxThrust = 0;
  double stage1MaxThrust = 0;
  double stage2MaxThrust = 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];
	stage1ExpansionRatio = x[4];
	stage2ExpansionRatio = x[10];

  if (stage == 2)
    {
      ExpansionRatio = x[10];
    }
  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;

    }


  stage2TankVolume = x[7] * (1 + (PayLoad / 1000));	//is as good as massratio ,but excludes otherdryweights
  stage2FuelMass =
    stage2TankVolume * FuelDensity * (1 - PressureTankFraction);
  for (stage = 1; stage <= NStages; stage++)
    {
      if (stage == 2)
	{
	stagingPenalty = 0; 
	  ExpansionRatio = stage2ExpansionRatio;  
	  //ExpansionRatio = x[4];  
	  double FlatVelocity = x[5];
#ifdef PF			//PRESSUREFED //defined above
	  double TankPressE = x[12];
#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

	  TankRadius = RadiusFrom (stage2TankVolume, stage2AspectRatio);
      AlternateTankThickness = ((TankPress * 2 * TankRadius) / MaterialStrength);	// standard formula stress = Pressure*Diameter 
      ActualTankThickness =
	fmax (MinimalMaterialThickness, AlternateTankThickness);
	stage2TankThickness = ActualTankThickness;

	  TankCrossSection = Pi * (TankRadius * TankRadius);
	  CylinderHeight = 2 * TankRadius * stage2AspectRatio;
	 stage2CylinderHeight = CylinderHeight;
	stage2TankRadius = TankRadius;
	  TankCircumference = 2 * Pi * TankRadius;
//double TankCrossSection =  ((TankRadius*TankRadius)*Pi);
	  stageTankVolume = stage2TankVolume;
	  stageFuelMass = stage2FuelMass;
	  TankStructureMass = (CylinderHeight * ActualTankThickness * TankCircumference * TankDensity + 6 * ActualTankThickness * TankCrossSection * TankDensity);	//  cylindrical tanks with spherical tops and bottoms   and a common bulkhead
//printf("ActualTankThickness = %G \n",ActualTankThickness);
	  stage2TankMass = TankStructureMass;
#ifdef PF
	  //EngineMass = (TankMass *etweightratio); 
	  EngineMass = (TankStructureMass + stageFuelMass) * stage2etweightratio;
	  //EngineMass = (StartFuelMass*etweightratio); 
	  //EngineMass = etweightratio; 
	  //EngineMass = 30; 
//printf("PF enginemass = %G etweightratio = %g tankmass = %g\n",EngineMass, etweightratio, TankMass);
#else
	  EngineMass = TankStructureMass * stage2etweightratio;
#endif
	stage2EngineMass = EngineMass;
	  //TankMass =  (TankStructureMass +RequiresLiner?TankLinerMass:0);
	  CurrentVehicleMass = stage2VehicleMass;
#define stage2EngineCost  (fmax(EngineCostPerWeight*0.,(stage2EngineMass *EngineCostPerWeight)))
#define stage2TankCost TankStructureMass*MaterialPrice
	  DeviceHeight = stage1CylinderHeight + stage2CylinderHeight + 2*( stage1TankRadius + stage2TankRadius);
	}
      if (stage == 1)
	{
	stagingPenalty = PayLoad/10; //this is for staging overhead
	printableStagingPenalty = stagingPenalty;
	//stagingPenalty = 0; //this is for staging overhead
	//stagingPenalty = 10; //this is for staging overhead

	  ExpansionRatio = stage1ExpansionRatio;  
	  //DeviceHeight = HeightFrom (stage1TankVolume, stage1AspectRatio);
	  TankRadius = RadiusFrom (stage1TankVolume, stage1AspectRatio);
      AlternateTankThickness = ((TankPress * 2 * TankRadius) / MaterialStrength);	// standard formula stress = Pressure*Diameter 
      ActualTankThickness =
	fmax (MinimalMaterialThickness, AlternateTankThickness);
	stage1TankThickness = ActualTankThickness;

	  TankCrossSection = Pi * (TankRadius * TankRadius);
	  CylinderHeight = 2 * TankRadius * stage1AspectRatio;
	 stage1CylinderHeight = CylinderHeight;
	stage1TankRadius = TankRadius;
	  TankCircumference = 2 * Pi * TankRadius;
//double TankCrossSection =  ((TankRadius*TankRadius)*Pi);
	  stageTankVolume = stage1TankVolume;
	  stageFuelMass =
	    stageTankVolume * FuelDensity * (1 - PressureTankFraction);
	  stage1FuelMass = stageFuelMass;
	  TankStructureMass = (CylinderHeight * ActualTankThickness * TankCircumference * TankDensity + 6 * ActualTankThickness * TankCrossSection * TankDensity);	//  cylindrical tanks with spherical tops and bottoms   and a common bulkhead

#ifdef PF
	  //EngineMass = (TankMass *etweightratio); 
	  EngineMass = (TankStructureMass + stageFuelMass) * stage1etweightratio;
	  //EngineMass = (stageFuelMass*etweightratio); 
	  //EngineMass = etweightratio; 
	  //EngineMass = 30; 
//printf("PF enginemass = %G etweightratio = %g tankmass = %g\n",EngineMass, etweightratio, TankMass);
#else
	  EngineMass = TankStructureMass * stage1etweightratio;
#endif
	stage1EngineMass = EngineMass;
	  //TankMass =  (TankStructureMass +RequiresLiner?TankLinerMass:0);
	  stage1TankMass = TankStructureMass;
	  CurrentVehicleMass = TotalVehicleMass;
	}

/* initialization stuff that depends on material  */
/* these used to be macros and rere recalculated 1000s of times */
      StartFuelMass = stageFuelMass;
      //TankLinerMass =(CylinderHeight*TankLinerThickness*TankCircumference*TankLinerDensity+6*TankLinerThickness*TankCrossSection*TankLinerDensity);


      //      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;


      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 (Burnout && stage < 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;// zzz  this is fucked for stage2
	  VMass = CurrentVehicleMass;
	  //FPress=TankBottomPress;
	  Thrust = CalcThrust;
	  if(Thrust >MaxThrust){
		if(stage == 1){
			stage1MaxThrust = Thrust;
		}else {
			stage2MaxThrust = Thrust;
		}

	  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 MFR %G\n",

		 NSec, X, Y, Acc, stageFuelMass, Theta, VelocityX, VelocityY,
		 Thrust, Drag,MFR);
	    }
	  //Acc = ThrustAndMassToAcc ;
	  stageFuelMass -= MFR;
	      CurrentVehicleMass -= MFR;
	  if (/*!Burnout && */(stageFuelMass < MFR || (MFR == 0)))
	    {
	      OutOfFuel = 1;
	      Burnout = NSec;
	      BurnoutAltitude = Y;
	      BurnoutVelocity = Velocity;
		if(stage == 1){
			break;
		}
	    }
	  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;
	  }
	}
	MaxThrust = 0;
	OutOfFuel = 0;
//printf("foobar \n");
//printf("###### stage %d  loop \n", stage);
//printf("barfoo \n");
    }
  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  \n",
	      CylinderHeight, TankCircumference, FlatVelocity);
      printf ("Profit = %G \n", Profit);
      printf ("MadeOrbit = %d\n", MadeOrbit);
      printf ("UninsuredLaunchCost = %G \n", UninsuredLaunchCost);
      printf ("YvelocityBogusfactor = %G\n", Yvelocitybogusfactor);
      printf ("Y = %G", Y);
      printf ("aspectVogosityFactor =    %G\n", aspectBogosityFactor);
      printf ("OrbitAltitude = %G \n", OrbitAltitude);
      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);
      //printf ("CylinderHeight=  %Gm \n", CylinderHeight);
      printf ("stage1AspectRatio = %G \n", stage1AspectRatio);
      printf ("stage2AspectRatio = %G \n", stage2AspectRatio);
      printf ("stage1TankMass = %G \n", stage1TankMass);
      printf ("stage2TankMass = %G \n", stage2TankMass);
      printf ("stage1TankCost = %G \n", stage1TankCost);
      printf ("stage2TankCost = %G \n", stage2TankCost);
      printf ("stage1FuelMass = %G \n", stage1FuelMass);
      printf ("stage2FuelMass = %G \n", stage2FuelMass);
	printf("stage1CylinderHeight = %G\n",stage1CylinderHeight);
	printf("stage1TankRadius = %G\n",stage1TankRadius);
	printf("stage2CylinderHeight = %G\n",stage2CylinderHeight);
	printf("stage2TankRadius = %G \n",stage2TankRadius);
	printf("stage1EngineMass = %G \n",stage1EngineMass);
	printf("stage2EngineMass = %G \n\n",stage2EngineMass);
	printf("stage1EngineCost = %G \n",stage1EngineCost);
	printf("stage2EngineCost = %G \n",stage2EngineCost);
	printf("TotalEngineCost = %G \n\n",TotalEngineCost);

	printf("stage1MaxThrust = %G \n",stage1MaxThrust);
	printf("stage2MaxThrust = %G \n",stage2MaxThrust);
	printf("stage1TankThickness = %G\n",stage1TankThickness);
	printf("stage2TankThickness = %G\n",stage2TankThickness);
	printf("stage1ExpansionRatio = %g \n",stage1ExpansionRatio);
	printf("stage2ExpansionRatio = %g \n",stage2ExpansionRatio);
	printf("stagingPenalty = %G\n",printableStagingPenalty);

	printf("MaterialDensity = %G \n",MaterialDensity);
      printf ("Device Height %Gm %Gft+ engine and Cone  diameter = %Gm %Gft\n",
	 DeviceHeight,
	 DeviceHeight * 3.280840, stage1TankRadius * 2,
	 stage2TankRadius * 2 * 3.280840);
      printf ("\naspect height/diameter %G\n",
	      DeviceHeight / (stage1TankRadius * 2));
//printf("\ncalculated approx MaxThrust %G N  %G lb\n",MaxThrust,NewtonsToLbf(MaxThrust));
/*zzz      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 ("TotalFuelMass   %Gkg  ", TotalFuelMass);
      printf ("TankMass %Gkg ", TotalTankMass);
      printf ("EngineMass %Gkg ", TotalEngineMass);
/*zzz      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);




 /*zzz     printf
	("DryMass  %Gkg DryMass+PayLoad %Gkg  PayLoad %Gkg StartFuelMass %Gkg  \nVehicle_Loaded  %G\n",
	 VehicleDryMass, VehicleDryMass + PayLoad, PayLoad, StartFuelMass,
	 VehicleDryMass + PayLoad + StartFuelMass);*/
/*zzz      printf ("MassRatio  %G   DeltaV   %G  \n\n",
	      (VehicleDryMass + PayLoad + StartFuelMass) / (VehicleDryMass +
							    PayLoad), DeltaV);
*/
   /*zzz   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));
//zzz      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("VehicleCost  $%G    Un//insuredLaunchCost  $%G\n",VehicleCost,UninsuredLaunchCost);
      printf ("\nVehicleCost  $%G    Tank $%G Engine $%G Fuel  $%G\n",
	      VehicleCost, TotalTankCost, TotalEngineCost, TotalFuelCost);
      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));
      //zzzprintf (" 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 / (TotalTankMass + TotalEngineMass),
	      (TotalTankMass + TotalEngineMass) / 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;
}

// 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;
}

/*double
  merit(){

	return Merit;
}*/