#ifndef COOLING_METAL_HINCLUDED
#define COOLING_METAL_HINCLUDED

/*
 * Cooling code for cosmology simulations.
 * Originally written by James Wadsley and Sijing Shen, McMaster
 * University for GASOLINE.
 */

/* Global consts */
/*#if defined(COOLDEBUG)
#include "mdl.h"
#endif*/
#include "param.h"
#include "rpc/xdr.h"

#ifdef __cplusplus
extern "C" {
#endif

#include "stiff.h"

/* Constants */
#define CL_B_gm         (6.022e23*(938.7830/931.494))/*Avegadro's Number * Mass_Hydrogen/Energy_AMU */
#define CL_k_Boltzmann  1.38066e-16
#define CL_eV_erg       1.60219e-12
#define CL_eV_per_K     (CL_k_Boltzmann/CL_eV_erg)
/*
#define CL_RT_FLOAT      float
#define CL_RT_MIN        1e-38
#define CL_RT_MIN        FLT_MIN
*/

#define CL_RT_FLOAT      double
#define CL_RT_MIN        1e-100

/*
#define CL_RT_MIN        DBL_MIN
*/
/* 
 * Work around for Dec ev6 flawed
 * treatment of sub-normal numbers 
 */
#define CL_MAX_NEG_EXP_ARG  -500.

#define CL_NMAXBYTETABLE   56000
#define MU_METAL  17.6003
#define ZSOLAR    0.0130215

typedef struct CoolingParametersStruct {
  int    bIonNonEqm;
  int    nCoolingTable;
  int    bUV;
  int    bMetal;
  char   *CoolInFile; 
  int    bUVTableUsesTime;
  int    bDoIonOutput;
  int    bLowTCool;
  int    bSelfShield;
  int    bShieldHI; /* Set to true if dust shields HI;*/ 
  double dMassFracHelium;
  double dCoolingTmin;     
  double dCoolingTmax;
  double dClump;
  double dLymanWernerFrac; /* Fraction of Lyman Werner radiation that escapes birth cloud.  0.5 is a good value.*/
} COOLPARAM;

typedef struct CoolingParticleStruct {
  double f_HI,f_HeI,f_HeII;
  double f_H2;	/* Abundance of ions */
  double     dLymanWerner; /* Flux of Lyman Werner radiation at the gas particle */
} COOLPARTICLE;

typedef struct { 
  double e,Total;
  double HI,HII,HeI,HeII,HeIII;
  double H2; 
} PERBARYON;

typedef struct { 
  double   zTime;

  double   Rate_Phot_HI;
  double   Rate_Phot_HeI;
  double   Rate_Phot_HeII;
  double   Rate_Phot_H2_cosmo; /* Dissociating radiation from the cosmic background for H2*/ 

  double   Heat_Phot_HI;
  double   Heat_Phot_HeI;
  double   Heat_Phot_HeII;
  double   Heat_Phot_H2; 
} UVSPECTRUM;

typedef struct { 
  double   Rate_Phot_HI;
  double   Rate_Phot_HeI;
  double   Rate_Phot_HeII;
  double   Rate_Phot_H2_cosmo;  

  double   Heat_Phot_HI;
  double   Heat_Phot_HeI;
  double   Heat_Phot_HeII;
  double   Heat_Phot_H2;  
 
  double   Cool_Coll_HI;
  double   Cool_Coll_HeI;
  double   Cool_Coll_HeII;
  double   Cool_Diel_HeII;
  double   Cool_Coll_H2;  
 
  double   Cool_Comp;
  double   Tcmb;
  double   Cool_LowTFactor;

} RATES_NO_T;

typedef struct { 
  CL_RT_FLOAT   Rate_Coll_HI;
  CL_RT_FLOAT   Rate_Coll_HeI;
  CL_RT_FLOAT   Rate_Coll_HeII;
  CL_RT_FLOAT   Rate_Coll_e_H2;  
  CL_RT_FLOAT   Rate_Coll_HI_H2;  
  CL_RT_FLOAT   Rate_Coll_H2_H2;  
  CL_RT_FLOAT   Rate_Coll_Hm_e;           /*gas phase formation of H2 */
  CL_RT_FLOAT   Rate_Coll_HI_e;           /*--------------------*/
  CL_RT_FLOAT   Rate_Coll_HII_H2;          /*--------------------*/
  CL_RT_FLOAT   Rate_Coll_Hm_HII;        /*-------------------- */
  CL_RT_FLOAT   Rate_HI_e;          /*-------------------- */
  CL_RT_FLOAT   Rate_HI_Hm;          /*gas phase formation of H2 */
  CL_RT_FLOAT   Rate_Radr_HII;
  CL_RT_FLOAT   Rate_Radr_HeII;
  CL_RT_FLOAT   Rate_Radr_HeIII;
  CL_RT_FLOAT   Rate_Diel_HeII;
  CL_RT_FLOAT   Rate_Chtr_HeII;

  CL_RT_FLOAT   Cool_Brem_1;
  CL_RT_FLOAT   Cool_Brem_2;
  CL_RT_FLOAT   Cool_Radr_HII;
  CL_RT_FLOAT   Cool_Radr_HeII;
  CL_RT_FLOAT   Cool_Radr_HeIII;
  CL_RT_FLOAT   Cool_Line_HI;
  CL_RT_FLOAT   Cool_Line_HeI;
  CL_RT_FLOAT   Cool_Line_HeII;
  CL_RT_FLOAT   Cool_Line_H2_H;   
  CL_RT_FLOAT   Cool_Line_H2_H2; 
  CL_RT_FLOAT   Cool_Line_H2_He; 
  CL_RT_FLOAT   Cool_Line_H2_e; 
  CL_RT_FLOAT   Cool_Line_H2_HII;  
  CL_RT_FLOAT   Cool_LowT;
} RATES_T;

typedef struct clDerivsDataStruct clDerivsData;

/* Heating Cooling Context */

typedef struct CoolingPKDStruct { 
  double     z; /* Redshift */
  double     dTime;
 /* Rates independent of Temperature */ 
  RATES_NO_T  R;
 /* Table for Temperature dependent rates */ 
  int        nTable;
  double     TMin;
  double     TMax;
  double     TlnMin;
  double     TlnMax;
  double     rDeltaTln;
  RATES_T     *RT;
  
  int         bMetal; 
  int         nzMetalTable;
  int         nnHMetalTable;
  int         nTMetalTable;  
  double      MetalTMin; 
  double      MetalTMax; 
  double      MetalTlogMin;
  double      MetalTlogMax;
  double      rDeltaTlog;
  double      MetalnHMin; 
  double      MetalnHMax; 
  double      MetalnHlogMin; 
  double      MetalnHlogMax; 
  double      rDeltanHlog;
  double      MetalzMin; 
  double      MetalzMax;  
  double      rDeltaz;
  float       ***MetalCoolln;
  float       ***MetalHeatln;  
  double      *Rate_DustForm_H2; 
  
  int        nTableRead; /* number of Tables read from files */

  int        bUV;
  int        nUV;
  UVSPECTRUM *UV;
  int        bUVTableUsesTime;
  int        bUVTableLinear;
  int        bLowTCool;
  int        bSelfShield;
  
  int        bShieldHI;
  double     dClump; /* Subgrid clumping factor for determining rate of H2 formation on dust.  10 is a good value*/
  double     dLymanWernerFrac; /*  Set to true to determine age of star particle from mass compared to formation mass when calculating LW radiation.  Useful in running ICs which already have stars*/
  double     dGmPerCcUnit;
  double     dComovingGmPerCcUnit;
  double     dExpand; /*cosmological expansion factor*/
  double     dErgPerGmUnit;
  double     dSecUnit;
  double     dErgPerGmPerSecUnit;
  double     diErgPerGmUnit;
  double     dKpcUnit;
  double     dMsolUnit;
  double     dMassFracHelium;

/* Diagnostic */
  int       its;
#if defined(COOLDEBUG)
  /*  MDL        mdl; *//* For diag/debug outputs */
  /*struct particle *p;*/ /* particle pointer needed for SN feedback */
  int        iOrder;
#endif 
} COOL;

typedef struct {
  double   T, Tln;
  double   Coll_HI;
  double   Coll_HeI;
  double   Coll_HeII;
  double   Coll_e_H2;  
  double   Coll_HI_H2;  
  double   Coll_H2_H2; 
  double   Coll_Hm_e;           /*gas phase formation of H2 */
  double   Coll_Hm_HII;          /*------------------- */
  double   Coll_HI_e;           /*------------------- */
  double   Coll_HII_H2;          /*--------------------- */
  double   HI_e;          /*---------------------- */
  double   HI_Hm;          /*gas phase formation of H2 */
  double   Radr_HII;
  double   Radr_HeII;
  double   Diel_HeII;
  double   Chtr_HeII;
  double   Totr_HeII;
  double   Radr_HeIII; 
  double   Cool_Metal;
  double   Heat_Metal;

  double   Phot_HI;
  double   Phot_HeI;
  double   Phot_HeII;
  double   Phot_H2;  /*Photon dissociation of H2*/
  double   DustForm_H2; /* Formation of H2 on dust */
  double   CorreLength; /* The correlation length of subgrid turbulence, used when calculating shielding*/
  double   LymanWernerCode;
} RATE;

typedef struct {
  double compton;
  double bremHII;
  double bremHeII;
  double bremHeIII;
  double radrecHII;
  double radrecHeII;
  double radrecHeIII;
  double collionHI; 
  double collionHeI;
  double collionHeII;
  double collion_e_H2;  
  double collion_H_H2;  
  double collion_H2_H2;
  double collion_HII_H2;
  double dielrecHeII;
  double lineHI;
  double lineHeI;
  double lineHeII;
  double lineH2;  
  double lowT;
  double NetMetalCool; 
} COOL_ERGPERSPERGM;


struct clDerivsDataStruct {
  STIFF *IntegratorContext;
  COOL *cl;
  double rho,ExternalHeating,E,ZMetal,dLymanWerner, columnL;
/*  double Y_H, Y_He; */  /* will be needed -- also for temperature , Y_MetalIon, Y_eMetal */
  RATE Rate;
  PERBARYON Y;
  double     Y_H, Y_He, Y_eMax;
  double     Y_Total0, Y_Total1;
  double     dlnE;
  int        its;  /* Debug */
  int        bCool;
};

COOL *CoolInit( );
void CoolFinalize( COOL *cl );
clDerivsData *CoolDerivsInit(COOL *cl);
void CoolDerivsFinalize(clDerivsData *cld ) ;

void clInitConstants( COOL *cl, double dGMPerCcunit, double dComovingGmPerCcUnit,
		      double dErgPerGmUnit, double dSecUnit, double dKpcUnit, COOLPARAM CoolParam);
void clInitUV(COOL *cl, int nTableColumns, int nTableRows, double *dTableData );
void clInitRatesTable( COOL *cl, double TMin, double TMax, int nTable );
void clReadMetalTable(COOL *cl, COOLPARAM clParam);
void clRateMetalTable(COOL *cl, RATE *Rate, double T, double rho, double Y_H, double ZMetal); 
void clHHeTotal(COOL *cl, double ZMetal); 
void CoolInitRatesTable( COOL *cl, COOLPARAM CoolParam);

void clRatesTableError( COOL *cl );
void clRatesRedshift( COOL *cl, double z, double dTime );
double clHeatTotal ( COOL *cl, PERBARYON *Y, RATE *Rate, double rho, double ZMetal );
void clRates( COOL *cl, RATE *Rate, double T, double rho, double ZMetal, double columnL, double  Rate_Phot_H2_stellar);
double clCoolTotal( COOL *cl, PERBARYON *Y, RATE *Rate, double rho, double ZMetal );
COOL_ERGPERSPERGM  clTestCool ( COOL *cl, PERBARYON *Y, RATE *Rate, double rho );
void clPrintCool( COOL *cl, PERBARYON *Y, RATE *Rate, double rho );
void clPrintCoolFile( COOL *cl, PERBARYON *Y, RATE *Rate, double rho, double ZMetal, FILE *fp );

void clAbunds( COOL *cl, PERBARYON *Y, RATE *Rate, double rho, double ZMetal);
double clThermalEnergy( double Y_Total, double T );
double clTemperature( double Y_Total, double E );
double clSelfShield (double yH2, double h);
double clDustShield (double yHI, double yH2, double z, double h);
double clRateCollHI( double T );
double clRateCollHeI( double T );
double clRateCollHeII( double T );
double clRateColl_e_H2( double T );
double clRateColl_HI_H2( double T );
double clRateColl_H2_H2( double T );
double clRateColl_HII_H2(double T);
double clRateColl_Hm_e(double T);
double clRateColl_HI_e(double T);
double clRateColl_Hm_HII(double T);
double clRateHI_e(double T);
double clRateHI_Hm(double T);
double clRateRadrHII( double T );
double clRateRadrHeII( double T );
double clRateDielHeII( double T );
double clRateChtrHeII(double T);
double clRateRadrHeIII( double T );
double clCoolBrem1( double T );
double clCoolBrem2( double T );
double clCoolRadrHII( double T );
double clCoolRadrHeII( double T );
double clCoolRadrHeIII( double T );
double clCoolLineHI( double T );
double clCoolLineHeI( double T );
double clCoolLineHeII( double T );
double clCoolLineH2_table( double T );  
double clCoolLineH2_HI( double T );
double clCoolLineH2_H2( double T );
double clCoolLineH2_He( double T );
double clCoolLineH2_e( double T );
double clCoolLineH2_HII( double T );
double clCoolLowT( double T );
double clRateDustFormH2(double z, double clump); 
double clEdotInstant ( COOL *cl, PERBARYON *Y, RATE *Rate, double rho, 
		       double ZMetal, double *dEdotHeat, double *dEdotCool );
  void clIntegrateEnergy(COOL *cl, clDerivsData *clData, PERBARYON *Y, double *E, 
		       double ExternalHeating, double rho, double ZMetal, double dt, double columnL, double dLymanWerner  );
  void clIntegrateEnergyDEBUG(COOL *cl, clDerivsData *clData, PERBARYON *Y, double *E, 
		       double ExternalHeating, double rho, double ZMetal,  double dt );


void clDerivs(double x, const double *y, double *yheat,
	      double *ycool, void *Data) ;
  
void CoolAddParams( COOLPARAM *CoolParam, PRM );
void CoolLogParams( COOLPARAM *CoolParam, FILE *fp );
void CoolOutputArray( COOLPARAM *CoolParam, int, int *, char * );

#define COOL_ARRAY0_EXT  "HI"
double COOL_ARRAY0(COOL *cl, COOLPARTICLE *cp, double ZMetal);
double COOL_SET_ARRAY0(COOL *cl, COOLPARTICLE *cp,double ZMetal, double val);

#define COOL_ARRAY1_EXT  "HeI"
double COOL_ARRAY1(COOL *cl, COOLPARTICLE *cp, double ZMetal);
double COOL_SET_ARRAY1(COOL *cl, COOLPARTICLE *cp,double ZMetal, double val);

#define COOL_ARRAY2_EXT  "HeII"
double COOL_ARRAY2(COOL *cl, COOLPARTICLE *cp, double ZMetal);
double COOL_SET_ARRAY2(COOL *cl, COOLPARTICLE *cp,double ZMetal, double val);

#define COOL_ARRAY3_EXT  "H2"
double COOL_ARRAY3(COOL *cl, COOLPARTICLE *cp, double ZMetal);
double COOL_SET_ARRAY3(COOL *cl, COOLPARTICLE *cp,double ZMetal, double val);

double COOL_EDOT( COOL *cl_, COOLPARTICLE *cp_, double ECode_, double rhoCode_, double ZMetal_, double *posCode_, double columnL_ );
#define COOL_EDOT( cl_, cp_, ECode_, rhoCode_, ZMetal_, posCode_, columnL_) (CoolCodeWorkToErgPerGmPerSec( cl_, CoolEdotInstantCode( cl_, cp_, ECode_, rhoCode_, ZMetal_, posCode_ , columnL_)))

double COOL_COOLING( COOL *cl_, COOLPARTICLE *cp_, double ECode_, double rhoCode_, double ZMetal_, double *posCode_ , double columnL_);
#define COOL_COOLING( cl_, cp_, ECode_, rhoCode_, ZMetal_, posCode_, columnL_) (CoolCodeWorkToErgPerGmPerSec( cl_, CoolCoolingCode( cl_, cp_, ECode_, rhoCode_, ZMetal_, posCode_ , columnL_)))

double COOL_HEATING( COOL *cl_, COOLPARTICLE *cp_, double ECode_, double rhoCode_, double ZMetal_, double *posCode_, double columnL_ );
#define COOL_HEATING( cl_, cp_, ECode_, rhoCode_, ZMetal_, posCode_, columnL_) (CoolCodeWorkToErgPerGmPerSec( cl_, CoolHeatingCode( cl_, cp_, ECode_, rhoCode_, ZMetal_, posCode_ , columnL_))) 

void clSetAbundanceTotals(COOL *cl, double ZMetal, double *Y_H, double *Y_He, double *Y_eMAX);
void CoolPARTICLEtoPERBARYON(COOL *cl_, PERBARYON *Y, COOLPARTICLE *cp, double ZMetal);
void CoolPERBARYONtoPARTICLE(COOL *cl_, PERBARYON *Y, COOLPARTICLE *cp, double ZMetal);

double CoolLymanWerner(double dAge);

double CoolEnergyToTemperature( COOL *Cool, COOLPARTICLE *cp, double E, double ZMetal);
double CoolCodeEnergyToTemperature( COOL *Cool, COOLPARTICLE *cp, double E, double ZMetal);

/* Note: nod to cosmology (z parameter) unavoidable unless we want to access cosmo.[ch] from here */
void CoolSetTime( COOL *Cool, double dTime, double z );

double CoolCodeTimeToSeconds( COOL *Cool, double dCodeTime );

#define CoolCodeTimeToSeconds( Cool, dCodeTime ) ((Cool)->dSecUnit*(dCodeTime))

double CoolSecondsToCodeTime( COOL *Cool, double dTime );

#define CoolSecondsToCodeTime( Cool, dTime ) ((dTime)/(Cool)->dSecUnit)

double CoolCodeEnergyToErgPerGm( COOL *Cool, double dCodeEnergy );

#define CoolCodeEnergyToErgPerGm( Cool, dCodeEnergy ) ((Cool)->dErgPerGmUnit*(dCodeEnergy))

double CoolErgPerGmToCodeEnergy( COOL *Cool, double dEnergy );

#define CoolErgPerGmToCodeEnergy( Cool, dEnergy ) ((Cool)->diErgPerGmUnit*(dEnergy))

double CoolCodeWorkToErgPerGmPerSec( COOL *Cool, double dCodeWork );

#define CoolCodeWorkToErgPerGmPerSec( Cool, dCodeWork ) ((Cool)->dErgPerGmPerSecUnit*(dCodeWork))

double CoolErgPerGmPerSecToCodeWork( COOL *Cool, double dWork );

#define CoolErgPerGmPerSecToCodeWork( Cool, dWork ) ((dWork)/(Cool)->dErgPerGmPerSecUnit)

double CodeDensityToComovingGmPerCc( COOL *Cool, double dCodeDensity );

#define CodeDensityToComovingGmPerCc( Cool, dCodeDensity )  ((Cool)->dComovingGmPerCcUnit*(dCodeDensity))

void CoolIntegrateEnergy(COOL *cl, clDerivsData *cData, COOLPARTICLE *cp, double *E, 
			 double ExternalHeating, double rho, double ZMetal, double tStep, double columnL );

void CoolIntegrateEnergyCode(COOL *cl, clDerivsData *cData, COOLPARTICLE *cp, double *E, 
			     double ExternalHeating, double rho, double ZMetal, double *r, double tStep, double columnL );

void CoolDefaultParticleData( COOLPARTICLE *cp );

void CoolInitEnergyAndParticleData( COOL *cl, COOLPARTICLE *cp, double *E, double dDensity, double dTemp, double ZMetal);

/* Deprecated */
double CoolHeatingRate( COOL *cl, COOLPARTICLE *cp, double E, double dDensity, double ZMetal, double columnL);

double CoolEdotInstantCode(COOL *cl, COOLPARTICLE *cp, double ECode, 
			   double rhoCode, double ZMetal, double *posCode, double columnL );
double CoolCoolingCode(COOL *cl, COOLPARTICLE *cp, double ECode, 
		       double rhoCode, double ZMetal, double *posCode, double columnL );
double CoolHeatingCode(COOL *cl, COOLPARTICLE *cp, double ECode, 
		       double rhoCode, double ZMetal, double *posCode, double columnL );

void CoolCodePressureOnDensitySoundSpeed( COOL *cl, COOLPARTICLE *cp, double uPred, double fDensity, double gamma, double gammam1, double *PoverRho, double *c );

/* Note: gamma should be 5/3 for this to be consistent! */
#define CoolCodePressureOnDensitySoundSpeed( cl__, cp__, uPred__, fDensity__, gamma__, gammam1__, PoverRho__, c__ ) { \
  *(PoverRho__) = ((5./3.-1)*(uPred__)); \
  *(c__) = sqrt((5./3.)*(*(PoverRho__))); }

/*
double CoolCodePressureOnDensity( COOL *cl, COOLPARTICLE *cp, double uPred, double fDensity, double gammam1 );

#define CoolCodePressureOnDensity( cl, cp, uPred, fDensity, gammam1 ) ((gammam1)*(uPred))
*/

void CoolTableReadInfo( COOLPARAM *CoolParam, int cntTable, int *nTableColumns, char *suffix );

void CoolTableRead( COOL *Cool, int nData, void *vData);

#ifdef __cplusplus
}
#endif

#endif