/*
This code compute  several feature of the habitable zone (HZ) and of the 4 Gyr
continuously habitable zone (CHZ), following the analytical models given in:  

Valle G., Dell'Omodarme M., Prada Moroni P.G., Degl'Innocenti S. (2013).
Evolution of the habitable zone of low-mass stars. Detailed stellar models and
analytical relationships for different mass and chemical compositions. 


To compile the code (with gcc):
gcc hz.c -lm -o HZ

The input file named "HZ-input.dat" must contain in each row:
1) the mass of the host star in solar unit (range 0.70 - 1.10)
2) the metallicity of the star (range: 0.005 - 0.04).
3) the initial helium content of the star.
4) the equilibrium temperature of the outer boundary (range: 169 - 203 K).
5) the albedo of the planet. 

The output is printed on screen.
 */


#include<stdio.h>
#include<stdlib.h>
#include<math.h>

main()
{
  FILE *fp;
  
  double M, Z, Y, T, T169, A;
  double K1, K2;
  double M2, M3, Ml, Zl, Zl2, Zl3, T169_2;
  
  double dm, tm, W, Ri, Ro;
  
  // coefficients for models computations
  double m_dm[] = {
    -1.24104436991515, 1.58195125658449, 1.07909978390663,
    1.08315938773256, -1.22895806944626, 0.095151118249582,
    0.547044851745198, 0.482169979415731, 0.867295533485984,
    -1.0128827071175, 0.00514252604909133, -0.279922101391439,
    -0.835057221157111, -0.0285136797071667, -0.0232323601030248 };

  double m_tm[] = {
    -37.2982440392992, 109.973591757898, 0.0664872690838029,
    -4.76605028487226, -36.3864322384106, 1.50978047217856,
    77.5872809671625, 40.7217089425163, 98.7414578768803,
    -0.855155227052455, -0.187062890781382, -212.544459415195,
    35.5695054528837, -76.4845080183614, -3.19633622253578,
    0.276434437355178,  0.319646320330445, -0.196686138381493,
    2.13735092885565 };

  double mod_w[] = {
    0.678921684261947, -0.976372657220196, -2.64867113446985,
    0.413812163186086, -0.234716920493049, 1.17403600475692,
    0.0941597565128401, 0.337098749956035, -1.70810899207781,
    -0.00651374005107242, 3.08588231502117, 2.12270394847918,
    1.61705242780183, -2.24114400998328 };

  double mod_ri[] = {
    -1.1957124470928, -0.271228308368989, -0.903647820092755,
    4.89475134379831, -0.244237810274819, -5.20397786462455,
    2.77387380477697, 0.731797295451257, -1.30592874210727,
    -3.6519677240594, 20.9248522519673, -10.244470379521,
    1.43831475362848,  6.03353002476469, 6.69106285185334,
    -0.1097577868901, -5.4294565822101, 2.01485572687134,
    10.5994415118025 };

  double mod_ro[] = {
    1.27358667250932, -1.67316177073267, -1.74649327615352,
    0.200947047137099, -1.19578441065701, 0.0479394883903738,
    2.36500799134778, 0.333568404450957, 0.236593021504852,
    -0.835587016321661, -0.329152934600133, 7.59377689532411,
    0.301779996733309, 0.382656976320151, -0.114535833536909,
    0.0703547942121432, -0.126327413929569, -8.9430362480413,
    3.17658310523942, -0.0109338842518824, -0.257222721955351 };

  // read input file
  fp = fopen("HZ-input.dat", "r");
  if(fp == NULL)
    {
      printf("Missing input file HZ-input.dat\n");
      exit(-1);
    }

  // output header
  printf("   dm     tm      W     Ri       Ro\n");
  printf("  (AU)  (Gyr)   (AU)   (AU)     (AU)\n");

  // cycle over the input rows
  while(fscanf(fp,"%lf %lf %lf %lf %lf", &M, &Z, &Y, &T, &A) != EOF)
    {
      T169 = T/169;
      T169_2 = T169*T169;

      M2 = M*M;
      M3 = M2*M;
      Ml = log10(M);
      
      Zl = log10(Z);
      Zl2 = Zl*Zl;
      Zl3 = Zl2*Zl;
  
      K1 = log10(sqrt((1-A)/0.7)/T169_2);
      K2 = log10(sqrt((1-A)/0.7));
      
      // computation of dm
      dm = K1 + m_dm[0] + m_dm[1]*M + m_dm[2]*Zl + m_dm[3]*Y + m_dm[4]*M2 +
	m_dm[5]*Zl3 +  m_dm[6]*Zl2 + m_dm[7]*M3 + m_dm[8]*Y*M2 +
	m_dm[9]*M*Zl +  m_dm[10]*Y*Zl3 + m_dm[11]*M2*Zl2 +
	m_dm[12]*Y*M3 + m_dm[13]*Zl3*M3 + m_dm[14]*Y*Zl3*M3;
      
      // computation of tm
      tm = m_tm[0] + m_tm[1]*Ml + m_tm[2]*Zl2 + m_tm[3]*Y + m_tm[4]*T169_2+
	m_tm[5]*Zl + m_tm[6]*T169 + m_tm[7]/M + m_tm[8]*Ml*T169_2 + 
	m_tm[9]*Y*Zl + m_tm[10]*Zl2*T169_2 + m_tm[11]*Ml*T169 + 
	m_tm[12]*T169_2/M + m_tm[13]*T169/M + m_tm[14]*Ml*Y +
	m_tm[15]*Zl2*T169 + m_tm[16]*Y*T169_2 + m_tm[17]*Zl*T169_2+
	m_tm[18]*Ml*Y*T169_2;

      // computation of W
      W = K2 + mod_w[0] + mod_w[1]*M + mod_w[2]*T169 + mod_w[3]*Zl2+
	mod_w[4]*Y + mod_w[5]*M2 + mod_w[6]*Zl3 + mod_w[7]*Zl +
	mod_w[8]*M2*T169 + mod_w[9]*Zl3*T169 + mod_w[10]*M*T169 +
	mod_w[11]*M*Y + mod_w[12]*Y*T169 + mod_w[13]*M*Y*T169;
      
      // computation of Ri
      Ri = K2 + mod_ri[0] + mod_ri[1]*M + mod_ri[2]*Zl2 + mod_ri[3]*Y +
	mod_ri[4]*Zl3 + mod_ri[5]*M2 + mod_ri[6]*M3 + mod_ri[7]*Zl +
	mod_ri[8]*M2*Zl2 + mod_ri[9]*M*Y +  mod_ri[10]*Y*M2 +  mod_ri[11]*Y*M3+
	mod_ri[12]*Y*Zl3 +  mod_ri[13]*Y*Zl +  mod_ri[14]*Y*Zl2 + 
	mod_ri[15]*M3*Zl3 +  mod_ri[16]*M*Zl +  mod_ri[17]*Y*M2*Zl2+
	mod_ri[18]*M*Y*Zl;

      // computation of Ro
      Ro = K2 + mod_ro[0] + mod_ro[1]*M + mod_ro[2]*T169 + mod_ro[3]*Zl +
	mod_ro[4]*Y + mod_ro[5]*Zl3 + mod_ro[6]*M2 + mod_ro[7]*T169_2 +
	mod_ro[8]*Zl2 + mod_ro[9]*M3 +  mod_ro[10]*M*Zl +  mod_ro[11]*M*Y +
	mod_ro[12]*M*T169 +  mod_ro[13]*Y*Zl +  mod_ro[14]*M2*Zl2 + 
	mod_ro[15]*Y*T169 +  mod_ro[16]*M2*T169 +  mod_ro[17]*M2*Y +
	mod_ro[18]*Y*M3 + mod_ro[19]*M3*Zl3 + mod_ro[20]*M*Y*Zl;

      // output
      printf("%6.3f %6.2f %6.2f %6.2f %6.2f\n", pow(10,dm), pow(10,tm), 
	     pow(10,W), pow(10,Ri), pow(10, Ro));
    }
}