#include <math.h>
#include "ParallelGravity.h"
#include "feedback.h"
#include "starlifetime.h"

#define max(A,B) ((A) > (B) ? (A) : (B))
#define min(A,B) ((A) < (B) ? (A) : (B))

/* MAXIMUM MASS: 109 M_sol, lifetime= 3.26 Myr */
void Padova::CoefInit (double dMetals) const
{
  double logZ;
  double Z = dMetals;
  
  Z = min(Z, zmax);
  Z = max(Z, zmin);
  logZ = log10(Z);
  a0 = a00 + a01*logZ + a02*logZ*logZ;
  a1 = a10 + a11*logZ + a12*logZ*logZ;
  a2 = a20 + a21*logZ + a22*logZ*logZ;
}

double Padova::Lifetime(double dStarMass, double dMetals) const
{
  /* finds stellar lifetime in yr corresponding to stellar 
     mass in solar masses */
  double logStarMass, logLtime, Ltime;
  
  logStarMass = log10(dStarMass);
  CoefInit (dMetals);
  logLtime = a0 + a1 *logStarMass + a2*logStarMass*logStarMass;
  Ltime = pow(10.0, logLtime);
  
  return Ltime;
}

double Padova::StarMass(double dStarLtime, double dMetals) const
{
  /* finds stellar mass in solar masses corresponding to stellar 
     lifetime dStarTime in yr */
  double logStarMass, StarMass;
  double a, b, c;
  
  if(dStarLtime <= 0.0)		/* Time can be zero */
    return 1000;
  
  CoefInit (dMetals);
  c = a0;
  c -= log10(dStarLtime);
  b = a1;
  a = a2;
  if(b*b - 4*a*c < 0.0)		/* time is too small for fitting
				   formula */
    return 1000;
  
  logStarMass = (-b - sqrt(b*b - 4*a*c))/(2*a);
  StarMass = pow(10., logStarMass);
  return StarMass;
}


#if 0
#include "supernova.h"

int
main(int argc, char **argv)
{
    int nsamp;
    int i;
    double dlgm;
    double lgm;
    double Ntot;
    double Mtot;

    PDVAPARAM ppdva;
    SN sn;
    double dStarMass;
    double mass[30], time[30];
    PARTICLE *p;
    double dMetals = 0.0;
    
    assert(argc == 2);
    
    nsamp = atoi(argv[1]);
    dlgm = (2.0 + 1.0)/nsamp;
    snInitialize (&sn);
	//    Ntot = dMSCumNumber(&MSparam, 0.0);
	//    Mtot = dMSCumMass(&MSparam, 0.0);
    
	//    dSTLtimeMStar (ppdva, sn, dStarMass, p);
	//    dSTMStarLtime (ppdva, sn, dStarLtime, p);
    
    PadovaInitialize(&ppdva);
	//    mass = 8;
	//    printf ("first %g %g\n", mass, dSTLtimeMStar (ppdva, mass, p));

    printf("mass(M_sol) time[yr]\n");        
    for(i = 0; i < nsamp; i++) {
	
		//	lgm = -1 + i*dlgm;
		lgm = i*dlgm;

		mass[i] = pow(10.0, lgm);
        time[i] = dSTLtimeMStar (ppdva, mass[i], dMetals);
		printf("%g %g\n", mass[i], time[i]);        

		}

    printf("time[yr] mass(M_sol)\n");        
    for(i = 0; i < nsamp; i++) {

        mass[i] = dSTMStarLtime (ppdva, time[i], dMetals);
		printf("%g %g\n", mass[i], time[i]);        
		}
    
  return 0;
}
#endif