/***********************************************************************
* Adaptive Simulated Annealing (ASA)
* Lester Ingber <ingber@ingber.com>
* Copyright (c) 1993-1996 Lester Ingber.  All Rights Reserved.
* The LICENSE file must be included with ASA code.
***********************************************************************/
//extern int printf (FILE *, ...);
extern int _flsbuf (unsigned int, FILE *);



#include <stdio.h>

 /* $Id: user_cst.h,v 13.5 1996/08/16 10:06:27 ingber Exp ingber $ */

 /* user_cst.h for Adaptive Simulated Annealing */

 /* Note that this is a trimmed version of the ASA_TEST problem.
    A version of this cost_function with more documentation and hooks for
    various templates is in user.c. */

 /* If you use this file to define your cost_function (the default),
    insert the body of your cost function just above the line
    "#if ASA_TEST" below.  (The default of ASA_TEST is FALSE.)

    If you read in information via the asa_opt file (the default),
    define *parameter_dimension and
    parameter_lower_bound[.], parameter_upper_bound[.], parameter_int_real[.]
    for each parameter at the bottom of asa_opt.

    The minimum you need to do here is to use
    x[0], ..., x[*parameter_dimension-1]
    for your parameters and to return the value of your cost function.  */

#if HAVE_ANSI
double
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)
#else
double
cost_function (x,
	       parameter_lower_bound,
	       parameter_upper_bound,
	       cost_tangents,
	       cost_curvature,
	       parameter_dimension,
	       parameter_int_real,
	       cost_flag,
	       exit_code,
	       USER_OPTIONS)
     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;
#endif
{

  /* *** Insert the body of your cost function here. *** */

void checkMR(double,double,double);
#if ASA_TEST
// Aluminum foil tube wrapped in Dyneema 
// tankage for LOX + Propane 
//  with 20% drop to chamber pressure.  Dyneema/aluminum become
// stronger at cryogenic temperatures.  -- Jim Bowery
 #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

#define NewtonsPerLbf 4.448222  //Constant
#define NewtonsToLbf(x)  (x/NewtonsPerLbf) //Constant
#define PascalsPerPsi 6894.75 //Constant
#define PascalsToPsi(x) (x/PascalsPerPsi)
#define PsiToPascals(x) (x*PascalsPerPsi)
#define AluminumDensity (2702.0) //Constant
#define PET 0
#define PETDensity       (1370.)//Constant
#define PETStrength       (3.44379e8)
#define PETPrice	   (5.) // Constant guess
#define FUZZ .0001
#define Dyneema 1
#define DyneemaDensity 970 //  this IS right http://dns.toyobo.co.jp/e/seihin/dn/dyneema/karui.htm
#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
//  zzzzz note calculate thickness from dyneema strength instead of this hack
#define DyneemaPrice (40*2.2)	//$40/lb but this was for aramid, not dyneema//Parameter
#define DcWater (.97*.62)//Parameter
#define GEarth 9.80665//Constant
#define EarthRadius 6377500
#define FuelPress PsiToPascals(5000)//Parameter
#define ChamberPress (FuelPress*.80)//Parameter
#define NominalMaxISP  350 //Parameter
//#define CalcISP (323.0)//Parameter
//#define CalcISP (250.0)//Parameter
#define Fudge .95//CONSTANT
#define CalcISP ( CharisticISP* CF_ThrustCoefficient*Fudge)
#define CharisticISP 214.3 //from sutton for LOX&RP1 with a CF_Thrust_coefficient of 1.4  pp194, 60

#define KSpecificHeatRatio 1.3
//#define CF_ThrustCoefficient (chamberperssure,expandedpressure,exetrnalpressure,ksp,area_expansion_ratio) (sqrt((2*ksp*ksp/(ksp-1.))*((2./(k+1))^((k+1)/(k-1)))*(1-((expandedpressure/chamberpressure)^((k-1)/k)))) + ((expandedpressure-externalpressure)/chamberPressure)* area_expansion_ratio)
// ksp is the specific reat ratio as defined by KSpecificHeadRatio
//I won't bother to simplify this, since it is all constants, and
//the compiler can deal with completely!!! how nice for expressive
//purposes!
//#define CF_ThrustCoefficient  1.6 // for 4000 psi engine ad 15 psi goes up above 1.8
#define CF_ThrustCoefficient  1.7 // for 4000 psi engine ad 15 psi goes up above 1.8
#define ExpansionRatio 25 //Parameter corresponds to optimal at 15psi & 4000psi with k = 1.3 
//#define ExpansionRatio 40 //Parameter
#define MassFlow (OutOfFuel?0.:(min(MaxDesiredG*VehicleMass/CalcISP ,(min(MassFlowAllowed,(TankCrossSection/ExpansionRatio)*(ChamberPress/CharisticISP))))))// stop burning when out of fuel
#define CalcThrust (MassFlow*CalcISP*GEarth)

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


#define VehicleMass  (VehicleDryMass+FuelMass+PayLoad)
#define MeritVehicleCost (FuelCost+LaborFactor*(TankCost+EngineCost))
#define VehicleCost (FuelCost+(TankCost+EngineCost))
#define TankCost TankMass*MaterialPrice
#define EngineCostPerWeight 100.
#define EngineCost  (max(EngineCostPerWeight*20.,(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 InsuranceRate .4
#define min(x,y) ((x)<(y)?(x):(y))
#define max(x,y) ((x)>(y)?(x):(y))
//#define OrbitalVelocityAtAltitude 7700
//#define OrbitAltitude		  200000
#define OrbitalVelocityAtAltitude 7375
#define OrbitAltitude             482800  //300 miles
#define MadeOrbit  ((VelocityX>OrbitalVelocityAtAltitude)&&(Y>OrbitAltitude))//
#define Price (MadeOrbit?PayLoad*(700.0*2.2):0)//Parameter
#define Profit (Price-LaunchCost)
//#define Merit (MadeOrbit?(Profit/UninsuredLaunchCost-GBogosityFactor):(max(OrbitAltitude,Y) *min(OrbitalVelocityAtAltitude,VelocityX))-10e10)//
#define Merit (MadeOrbit?(Profit/UninsuredLaunchCost):(min(OrbitAltitude,max(1,Y)) *min(OrbitalVelocityAtAltitude,max(1,VelocityX))-10e10))//

//#define Merit (MadeOrbit?10e10-VehicleCost:(max(OrbitAltitude,Y)*min(OrbitalVelocityAtAltitude,VelocityX))-10e10)//
#define MinimalMaterialThickness .001//Parameter
#define	ActualTankThickness (max(MinimalMaterialThickness,AlternateTankThickness))
#define AlternateTankThickness  ((TankPress*2*TankRadius)/MaterialStrength) // standard formula stress = Pressure*Diameter
#define Pi (3.141592654)//Constant  !! it sure is!!
#define TankRadius (TankCircumference/(2.0*Pi))
#define	TankCrossSection ((TankRadius*TankRadius)*Pi)
#define TankVolume (TankHeight*TankCrossSection+(4/3)*Pi*(TankRadius*TankRadius*TankRadius)*2)
#define	TankDensity  MaterialDensity
#define MaxDesiredG 7
#define MaxDesiredAcc (MaxDesiredG*GEarth)
#define MaxAllowedAcc (100*GEarth) // in gs  this is actually about ok lets see what effect it has on the outcome, it's not actually a max its more of adesign spec and the multiplier is to compensate for execeeding it.//Parameter
#define MinDesiredAcc (2*GEarth)
//#define TankMass (TankStructureMass +TankLinerMass)//Dyneema
#define TankMass (TankStructureMass +RequiresLiner?TankLinerMass:0)
#define	TankStructureMass (TankHeight*ActualTankThickness*TankCircumference*TankDensity+8*ActualTankThickness*TankCrossSection*TankDensity+TankCircumference*2*TankRadius*MinimalMaterialThickness*TankDensity)//  cylindrical tanks with spherical tops and bottoms  with a cylinder connecting the two tanks     zzz add 1mm of aluminum to the insides of the tanks!
#define	TankLinerMass (TankHeight*TankLinerThickness*TankCircumference*TankLinerDensity+8*TankLinerThickness*TankCrossSection*TankLinerDensity)// 
//#define TankLinerMass 0
#define TankLinerDensity AluminumDensity
#define TankLinerThickness .001
#define	EngineThickness (AvgTankThickness*2.0)// just bogus
#define	EngineLength (4.0)//CONSTANT likewise bogus
#define	EngineDensity (2*AluminumDensity)
//#define	EngineMass (TankCircumference*EngineThickness*EngineLength*EngineDensity)
#define  ThrustWeightRatioOfEngine (100.) // 100 high normal range we might be 300 or 1000//Parameter
#define MinEngineMass 0.
#define NominalEngineMass (MaxThrust/(ThrustWeightRatioOfEngine*GEarth))
#define EngineMass (max(MinEngineMass,NominalEngineMass))
#define	VehicleDryMass (TankMass+EngineMass)
#define	FuelDensity (1100.0)//CONSTANT
#define CoefficientOfDrag 0.1//CONSTANT
#define Drag (CoefficientOfDrag/2.0)*AirDensity*Velocity*Velocity*TankCrossSection // what units is Drag ? zzz  newtons??
#define AirDensity (airdens[min(200,(int)(Y/1000.0))])//Cnstant
//#define FlightPath(x)  (max(0.,(AFlight +BFlight *(x) + CFlight*(x)*(x)+DFlight*(x)*(x)*(x))))
//#define FlightPath(x)  ((AFlight +BFlight *(x) + CFlight*(x)*(x)+DFlight*(x)*(x)*(x)))
#define FlightPath(x)  (FlatVelocity)
  static LONG_INT funevals = 0;

  double TankHeight=x[0];
  double TankCircumference=x[1];
  double MassFlowAllowed=x[2];
  double FlatVelocity=x[3];

  int MaterialType = x[4];//0 PET 1 Dyneema + liner
double PayLoad = 20 ;  // this is for the experimental rocket  //Parameter
//  double AFlight=x[4];
//  double BFlight=x[5];
//  double CFlight=x[6];
//  double DFlight=x[6];
 double MaxThrust=NominalMaxISP*MassFlowAllowed*GEarth;
//  double ApScale=x[6];
//  double PsiHigh=x[7];
//  double PsiLow=x[8];

/* Begin meat of cost fuction. */

// Initial conditions
//  double TankHeight = TankCircumference if TankHeight<TankCircumference;
  double Acc = GEarth;
  double VelocityX=0,VelocityY=0,Velocity=0,FPress=0;
  double Y = 0;
  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 MaxISP = 0;
// double MaxISPPressure = 0;
  double MaxY = 0, MaxV = 0;
  double airdens[300];
  double FuelMass = TankVolume*FuelDensity;
  double StartFuelMass=FuelMass;
static double MaxMerit=-10e18;
static int testNumber = 0;
static int printedTestNumber = 0;	
  int NSec = 0;
  double GBogosityFactor = 0;
  double TakeOffThrust = 0;
  //double FractionOfOrbital = 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;
 static int debugPrint = 0;
airdens[0]=1.2250;
airdens[1]=1.111;
airdens[2]=1.0666;
airdens[3]=.909;
airdens[4]=0.819;
airdens[5]=0.736;
airdens[6]=0.05;
airdens[7]=0.05;
airdens[8]=0.05;
airdens[9]=0.05;
airdens[10]=0.05;
airdens[11]=0.05;
airdens[12]=0.05;
airdens[13]=0.05;
airdens[14]=0.05;
airdens[15]=0.05;
airdens[16]=0.05;
airdens[17]=0.05;
airdens[18]=0.05;
airdens[19]=0.05;
airdens[20]=0.005;
airdens[21]=0.005;
airdens[22]=0.005;
airdens[23]=0.005;
airdens[24]=0.005;
airdens[25]=0.005;
airdens[26]=0.005;
airdens[27]=0.005;
airdens[28]=0.005;
airdens[29]=0.005;
airdens[30]=0.005;
airdens[31]=0.005;
airdens[32]=0.005;
airdens[33]=0.005;
airdens[34]=0.0005;
airdens[35]=0.0005;
airdens[36]=0.0005;
airdens[37]=0.0005;
airdens[38]=0.0005;
airdens[39]=0.0005;
airdens[40]=0.0005;
airdens[41]=0.0005;
airdens[42]=0.0005;
airdens[43]=0.0005;
airdens[44]=0.0005;
airdens[45]=0.0005;
airdens[46]=0.0005;
airdens[47]=0.0005;
airdens[48]=0.0005;
airdens[49]=0.0005;
airdens[50]=0.0005;
airdens[51]=0.0005;
airdens[52]=0.0005;
airdens[53]=0.0005;
airdens[54]=0.0005;
airdens[55]=0.0005;
airdens[56]=5e-05;
airdens[57]=5e-05;
airdens[58]=5e-05;
airdens[59]=5e-05;
airdens[60]=5e-05;
airdens[61]=5e-05;
airdens[62]=5e-05;
airdens[63]=5e-05;
airdens[64]=5e-05;
airdens[65]=5e-05;
airdens[66]=5e-05;
airdens[67]=5e-05;
airdens[68]=5e-05;
airdens[69]=5e-05;
airdens[70]=5e-06;
airdens[71]=5e-06;
airdens[72]=5e-06;
airdens[73]=5e-06;
airdens[74]=5e-06;
airdens[75]=5e-06;
airdens[76]=5e-06;
airdens[77]=5e-06;
airdens[78]=5e-06;
airdens[79]=5e-06;
airdens[80]=5e-07;
airdens[81]=5e-07;
airdens[82]=5e-07;
airdens[83]=5e-07;
airdens[84]=5e-07;
airdens[85]=5e-07;
airdens[86]=5e-07;
airdens[87]=5e-07;
airdens[88]=5e-07;
airdens[89]=5e-07;
airdens[90]=5e-07;
airdens[91]=5e-07;
airdens[92]=5e-07;
airdens[93]=5e-07;
airdens[94]=5e-07;
airdens[95]=5e-07;
airdens[96]=5e-07;
airdens[97]=5e-07;
airdens[98]=5e-07;
airdens[99]=5e-07;
airdens[100]=5e-08;
airdens[101]=5e-08;
airdens[102]=5e-08;
airdens[103]=5e-08;
airdens[104]=5e-08;
airdens[105]=5e-08;
airdens[106]=5e-08;
airdens[107]=5e-08;
airdens[108]=5e-08;
airdens[109]=5e-08;
airdens[110]=5e-09;
airdens[111]=5e-09;
airdens[112]=5e-09;
airdens[113]=5e-09;
airdens[114]=5e-09;
airdens[115]=5e-09;
airdens[116]=5e-09;
airdens[117]=5e-09;
airdens[118]=5e-09;
airdens[119]=5e-09;
airdens[120]=5e-09;
airdens[121]=5e-09;
airdens[122]=5e-09;
airdens[123]=5e-09;
airdens[124]=5e-09;
airdens[125]=5e-09;
airdens[126]=5e-09;
airdens[127]=5e-09;
airdens[128]=5e-09;
airdens[129]=5e-09;
airdens[130]=5e-10;
airdens[131]=5e-10;
airdens[132]=5e-10;
airdens[133]=5e-10;
airdens[134]=5e-10;
airdens[135]=5e-10;
airdens[136]=5e-10;
airdens[137]=5e-10;
airdens[138]=5e-10;
airdens[139]=5e-10;
airdens[140]=5e-10;
airdens[141]=5e-10;
airdens[142]=5e-10;
airdens[143]=5e-10;
airdens[144]=5e-10;
airdens[145]=5e-10;
airdens[146]=5e-10;
airdens[147]=5e-10;
airdens[148]=5e-10;
airdens[149]=5e-10;
airdens[150]=5e-10;
airdens[151]=5e-10;
airdens[152]=5e-10;
airdens[153]=5e-10;
airdens[154]=5e-10;
airdens[155]=5e-10;
airdens[156]=5e-10;
airdens[157]=5e-10;
airdens[158]=5e-10;
airdens[159]=5e-10;
airdens[160]=5e-10;
airdens[161]=5e-10;
airdens[162]=5e-10;
airdens[163]=5e-10;
airdens[164]=5e-10;
airdens[165]=5e-10;
airdens[166]=5e-10;
airdens[167]=5e-10;
airdens[168]=5e-10;
airdens[169]=5e-10;
airdens[170]=5e-11;
airdens[171]=5e-11;
airdens[172]=5e-11;
airdens[173]=5e-11;
airdens[174]=5e-11;
airdens[175]=5e-11;
airdens[176]=5e-11;
airdens[177]=5e-11;
airdens[178]=5e-11;
airdens[179]=5e-11;
airdens[180]=5e-11;
airdens[181]=5e-11;
airdens[182]=5e-11;
airdens[183]=5e-11;
airdens[184]=5e-11;
airdens[185]=5e-11;
airdens[186]=5e-11;
airdens[187]=5e-11;
airdens[188]=5e-11;
airdens[189]=5e-11;
airdens[190]=5e-11;
airdens[191]=5e-11;
airdens[192]=5e-11;
airdens[193]=5e-11;
airdens[194]=5e-11;
airdens[195]=5e-11;
airdens[196]=5e-11;
airdens[197]=5e-11;
airdens[198]=5e-11;
airdens[199]=5e-11;
airdens[200]=5e-11;
airdens[201]=5e-11;
airdens[202]=5e-11;
airdens[203]=5e-11;
airdens[204]=5e-11;
airdens[205]=5e-11;
airdens[206]=5e-11;
airdens[207]=5e-11;
airdens[208]=5e-11;
airdens[209]=5e-11;
airdens[210]=5e-11;
airdens[211]=5e-11;
airdens[212]=5e-11;
airdens[213]=5e-11;
airdens[214]=5e-11;
airdens[215]=5e-11;
airdens[216]=5e-11;
airdens[217]=5e-11;
airdens[218]=5e-11;
airdens[219]=5e-11;
airdens[220]=5e-11;
airdens[221]=5e-11;
airdens[222]=5e-11;
airdens[223]=5e-11;
airdens[224]=5e-11;
airdens[225]=5e-11;
airdens[226]=5e-11;
airdens[227]=5e-11;
airdens[228]=5e-11;
airdens[229]=5e-11;
airdens[230]=5e-11;
airdens[231]=5e-11;
airdens[232]=5e-11;
airdens[233]=5e-11;
airdens[234]=5e-11;
airdens[235]=5e-11;
airdens[236]=5e-11;
airdens[237]=5e-11;
airdens[238]=5e-11;
airdens[239]=5e-11;
	switch (MaterialType){
              case PET: //PET
                  MaterialPrice =PETPrice;	
                  MaterialStrength = PETStrength;
                  MaterialDensity = PETDensity;
                  MaterialName = "PET";
		  RequiresLiner = FALSE;
                  break;
               case Dyneema://Dyneema
                  MaterialPrice =DyneemaPrice;	
                  MaterialStrength = DyneemaStrength;
                  MaterialDensity = DyneemaDensity;
                  MaterialName = "Dyneema";
		  RequiresLiner = TRUE;
                  break;

               

}
	while(VelocityY>=-1.0 && Y>-1.){// we let VelocityY go negative because we havent tested  centrifugal force yet NOte add that  
		//if(MadeOrbit) {break;}
		if(VelocityX <-FUZZ){break;}
		if(Theta>(Pi/2)){break;}
		if(NSec >1000){break;}
		//if(Acc > MaxAllowedAcc) {break;}
		Acc = ThrustAndMassToAcc ;
		GBogosityFactor+=((( max(0,Acc -MaxDesiredAcc)+max(0,MinDesiredAcc -Acc))/GEarth)*PayLoad);  // integrate out of range acceleration and normalize wrt rocket size by factoring in payload
		VelocityX += cos(Theta)*Acc ;
		VelocityY += sin(Theta)*Acc - GEarth + VelocityX*VelocityX/(EarthRadius+Y);
		//VelocityY += sin(Theta)*Acc - GEarth(EarthRadius*EarthRadius/((EarthRadius + Y)*(EarthRadius +Y)) + VelocityX*VelocityX/(EarthRadius+Y);
		Velocity = sqrt(VelocityX*VelocityX+VelocityY*VelocityY);
		Y += VelocityY;
		X += VelocityX;
		DeltaV += Acc;
		MFR=MassFlow;
		VMass=VehicleMass;
		FPress=TankBottomPress;
		ISP=CalcISP;
		if(ISP >MaxISP){
			 MaxISP = ISP;

		}
		THRUSTSUM += ISP*MFR*GEarth;
		Thrust=CalcThrust;
		if(NSec == 1){TakeOffThrust = Thrust;}

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





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


//printf(" sec = %d X %G Y %G Theta %G VelocityX %G VelocityY %G  Velocity %G\n",NSec,X,Y,Theta,VelocityX,VelocityY,Velocity);
//NSec 	 X Y Acc FuelRemaining Theta VelocityX VelocityY
if(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(FuelMass <MFR) { OutOfFuel = 1;Burnout = NSec;}
		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;
	}
    *cost_flag = TRUE;
//  if(Y<OrbitAltitude||VelocityX<OrbitalVelocityAtAltitude) *cost_flag = FALSE;
/*if(PascalsToPsi(MaxFPress)>=2000.0) { 

	*cost_flag = FALSE;
	Y = X  = VelocityX = VelocityY = 0;
	
}
*/
testNumber +=1;
if(debugPrint){
	printf("TankHeight %G TankCircumference %G MassFlowAllowed %G FlatVelocity %G  Merit %G\n",TankHeight,TankCircumference,MassFlowAllowed,FlatVelocity,Merit); 
}
if(MaxMerit<Merit){ 
printedTestNumber +=1;
MaxMerit=Merit;
putchar(0x1b);
putchar(0x5b);
putchar(0x3b);
putchar(0x48);
putchar(0x1b);
putchar(0x5b);
putchar(0x32);
putchar(0x4a);

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

printf("   Flight path parameters \n");
printf("      FlatVelocity=  %G\n",FlatVelocity);
//printf("      AFlight %G",AFlight);
//printf("      BFlight %G",BFlight);
//printf("      CFlight %G",CFlight);
//printf("      DFlight %G",DFlight);
printf("\n");
printf("VelocityX =  %Gm/s\n Altitude=  %Gm  VelocityY %Gm/sec  Velocity %Gm/s\n",VelocityX,Y,VelocityY,Velocity);
printf("ISP=  %G   MaxAcc=  %Gg MinAcc %Gg   MaxFPress=  %Gpsi\n",THRUSTSUM/(StartFuelMass*GEarth),MaxAcc/GEarth,MinAcc/GEarth,PascalsToPsi(MaxFPress));
printf("MinAccAltitude %gm MinAccVelocity %Gm/sec \n",MinAccAltitude,MinAccVelocity);
printf("Burnout at %d seconds   Time T = %d sec\n",Burnout,NSec);

printf("\nA few constants\n");
printf("engine thrust/weight ratio %G\n",ThrustWeightRatioOfEngine);
printf("  PayLoad=  %Gkg\n",PayLoad);
printf("\ncost assumptions:\n  LaborFactor %G set to 1 for experimental \n  FuelUnitCost $%G/m^3\n  EnginePrice $%G/kg\n  %sPrice $%G/kg\n",LaborFactor,FuelUnitCost,EngineCostPerWeight,MaterialName,MaterialPrice);
printf("Thrust  %Gn %Glb   Thrust/Weight   %G \n",TakeOffThrust,NewtonsToLbf(TakeOffThrust	),TakeOffThrust/(GEarth*(VehicleDryMass+PayLoad+StartFuelMass)));

printf("\nStartFuelMass   %Gkg  \n",StartFuelMass);

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


//printf("FuelDensity   %G FuelMass   %G TankVolume* FuelDensity   %Gkg\n",FuelDensity, FuelMass, TankVolume * FuelDensity);
//printf("MaxISP = %G  MaxISPPressure = %G psi\n",MaxISP,PascalsToPsi(MaxISPPressure));
printf("MaxISP = %G  \n",MaxISP);


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


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

	checkMR( THRUSTSUM/(StartFuelMass*GEarth),(VehicleDryMass+PayLoad+StartFuelMass)/(VehicleDryMass+PayLoad),DeltaV);
//printf("VehicleCost  $%G    Un//insuredLaunchCost  $%G\n",VehicleCost,UninsuredLaunchCost);
printf("\nVehicleCost  $%G    TankCost $%G EngineCost $%G FuelCost  $%G\n",VehicleCost, TankCost, EngineCost, FuelCost);
printf("  MeritVehicleCost  with LaborFactor for the anealer %G\n",MeritVehicleCost);
 printf("MaxDrag=%G at Alt=%G\n",MaxDrag,MaxDragAlt);
printf("GBogosityFactor = %G\n",GBogosityFactor);
printf("Profit= %G\n\n",Profit);

printf("MadeOrbit  %s   Merit %G\n\n",MadeOrbit?"Yup":"NO GODDAMNIT!",Merit);
printf("\n");

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

  return (-Merit);
#endif /* ASA_TEST */
}

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


}