#include <math.h>
#include "ParallelGravity.h"
#include "starlifetime.h"
#include "supernova.h"
#include "romberg.h"

/* Calculation of number and mass of stars that go SN Type Ia in a
   given stellar mass range.  */

#define EPSSNIA 1e-6

/// class to hold information about the secondary
class MSBinary 
{
public:
    double dMass2;  // mass of secondary
    double dGamma_2nd;  // power law index of secondary distribution
    IMF *imf;
    };
    
/// Integrand for secondary mass function: probability of a binary
/// of mass dMassB having a secondary of MSBinary->dMass2
double dMSIMFSecInt(const MSBinary *msb, double dMassB)
{
    // A factor of 1/(dMassB*log(10)) is from the definition of IMF as
    // number per log10(M).  Another factor of 1/dMassB is to
    // normalized the integral over mass ratios.
    return pow(msb->dMass2/dMassB, msb->dGamma_2nd)
	*msb->imf->returnimf(dMassB)/dMassB/dMassB/log(10.0);
    }

/** IMF of secondary in binary system that goes SN Ia
 * The distribution of secondary mass ratios is assumed to be a power
 * law, mu^dGamma as in Matteucci & Greggio (1986)
 */
double dMSIMFSec(const SN *sn, double dMass2)
{
    double dIMFSec;
    double dMass2_2, Msup, Minf;
    const double dGamma_2nd = 2.0;
    const double dNorm_2nd = pow(2.0, 1 + dGamma_2nd)*(1 + dGamma_2nd);  // Normalization
									 // of 2ndary 

    MSBinary msb = {dMass2, dGamma_2nd, sn->imf};

    Msup = dMass2 + sn->dMSNIImin;  // Mass where primary would have gone supernovaII
    // Msup = 16.0;  // Note that the double integral can be tested by
    // uncommenting the above line, setting dFracBinSNIa = 1.0, and
    // comparing NSNIa with the number of stars in the 3.0-16.0 mass interval.
    dMass2_2 = 2.*dMass2;  // Minimum mass of binary
    Minf = (dMass2_2 > sn->dMBmin)?dMass2_2:sn->dMBmin;
    dIMFSec = dRombergO(&msb, (double (*)(const void *, double))dMSIMFSecInt, Minf, Msup, EPSSNIA);
    dIMFSec *= sn->dFracBinSNIa * dNorm_2nd;
    return dIMFSec;
    }

/** calculate number of SN Type Ia a la Raiteri, Villata, Navarro, A&A
  * 315, 105, 1996) Returns number of SN Type Ia that occur during
  * timestep in which dMassT1 and dMassT2 are masses of stars that end
  * their lives at the end and beginning of timestep, respectively
  */
double SN::NSNIa (double dMassT1, double dMassT2) const
{
    CkAssert (dMassT1 < dMassT2);
    // Exclude primaries that go SNII
    if (dMassT1 >= 0.5*dMBmax)
	return 0.0;
    if(dMassT2 > 0.5*dMBmax)
	dMassT2 = 0.5*dMBmax;
    return dRombergO(this, (double (*)(const void *, double))dMSIMFSec, dMassT1, dMassT2, EPSSNIA);
    }

#ifdef SNIA_TST
#include "testsnia.decl.h"
#define max(A,B) ((A) > (B) ? (A) : (B))
#define min(A,B) ((A) < (B) ? (A) : (B))
#include "feedback.h"

int
main(int argc, char **argv)
{
    Chabrier chimf; // Chabrier has power law index of -1.3 which
		    // matches s = -2.35 in Greggio & Renzini, 1983
    // Kroupa93 chimf; // Kroupa is what RVN use.
    // Kroupa01 chimf;
    // MillerScalo chimf;
    SN sn;
    sn.imf = &chimf;
    Padova pdva;
    double z = 0.02;		// metalicity: use value to compare
				// with Greggio & Renzini, 1983

    double nsamp = 100;
    double tfac = log(14.0e9/1e6)/nsamp;  /// equal log interavals from 1Myr
				    /// to 14Gyr
    
    // sn.dFracBinSNIa = 1.00;
    printf("# Total Ia, II supernova: %g %g, SNIa mass range: %g\n",
	   sn.NSNIa(0.8, 8.0)/chimf.CumMass(0.0),
	   (chimf.CumNumber(8.0) - chimf.CumNumber(40.0))/chimf.CumMass(0.0),
	   (chimf.CumNumber(3.0) - chimf.CumNumber(16.0))/chimf.CumMass(0.0));
    {
	// One solar mass formed at t = 0, with metallicity 0.02
	SFEvent sfEvent(1.0, 0.0, 0.02, .005, .005);
	FBEffects fbEffectsII;
	FBEffects fbEffectsIa;
	double t, deltat;
	double dMOxII = 0.0;
	double dMFeII = 0.0;
	double dMOxIa = 0.0;
	double dMFeIa = 0.0;
	
	deltat = 1e6;
	for (int i = 0; i < 10000; i++) {
	    t = i*deltat;
	    sn.CalcSNIIFeedback(&sfEvent, t, deltat, &fbEffectsII);
	    dMOxII += fbEffectsII.dMassLoss*fbEffectsII.dMOxygen;
	    dMFeII += fbEffectsII.dMassLoss*fbEffectsII.dMIron;
	    sn.CalcSNIaFeedback(&sfEvent, t, deltat, &fbEffectsIa);
	    dMOxIa += fbEffectsIa.dMassLoss*fbEffectsIa.dMOxygen;
	    dMFeIa += fbEffectsIa.dMassLoss*fbEffectsIa.dMIron;
	    }
	printf("# Total II Ox: %g, Fe: %g\n", dMOxII, dMFeII);
	printf("# Total Ia Ox: %g, Fe: %g\n", dMOxIa, dMFeIa);
	}
    
    for(int i = 0; i < nsamp; i++) {
	double t = 1e6*exp(i*tfac);
	double deltat = 0.01*t;  /// interval to determine rate
	double dMaxMass = pdva.StarMass(t, z); // stars dying at start
	double dMinMass = pdva.StarMass(t+deltat, z); // stars dying
						      // at end
	double dMStarMinII = max (8.0, dMinMass); // range of SNII
	double dMStarMaxII = min (40.0, dMaxMass);
	double dCumNMinII = chimf.CumNumber(dMStarMinII); 
	double dCumNMaxII = chimf.CumNumber(dMStarMaxII);
	double dMtot = chimf.CumMass(0.0);
	double dNSNII;
	// One solar mass formed at t = 0, with metallicity 0.02
	SFEvent sfEvent(1.0, 0.0, 0.02, .005, .005);
	FBEffects fbEffectsII;
	FBEffects fbEffectsIa;
	
	if(dMaxMass > 8.0 && dMinMass < 40.0) {
	    dNSNII = (dCumNMinII - dCumNMaxII)/dMtot/deltat;
	    }
	else dNSNII = 0.0;
	
	sn.CalcSNIIFeedback(&sfEvent, t, deltat, &fbEffectsII);
	sn.CalcSNIaFeedback(&sfEvent, t, deltat, &fbEffectsIa);
	
	if (dMaxMass > dMinMass)
	    printf ("%g %g %g %g %g %g %g %g %g %g %g %g %g\n", t,
		    sn.NSNIa(dMinMass, dMaxMass)/deltat/dMtot, dNSNII,
		    fbEffectsIa.dEnergy, fbEffectsII.dEnergy,
		    fbEffectsIa.dMassLoss, fbEffectsII.dMassLoss,
		    fbEffectsIa.dMetals, fbEffectsII.dMetals,
		    fbEffectsIa.dMIron, fbEffectsII.dMIron,
		    fbEffectsIa.dMOxygen, fbEffectsII.dMOxygen
		    );
	}
    }
#include "testsnia.def.h"
#endif