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

/* Command line tool to calculate HZ limits (Kopparapu et al. 2013).
 * 
 * Rory Barnes, 19 Mar 2013
 *
 * Input a star's luminosity and effective temperature to obtain the limits.
 * To compile: gcc -o hzlims hz2013.c -lm.
 * Usage: hzlims -l <luminosity> -t <teff>
 */

void HabitableZone(double l,double teff,double lim[]) {
  int i;
  double tstar,seff[6];
  double seffsun[6],a[6],b[6],c[6],d[6];

  seffsun[0] = 1.7763;
  seffsun[1] = 1.0385;
  seffsun[2] = 1.0146;
  seffsun[3] = 0.3507;
  seffsun[4] = 0.2946;
  seffsun[5] = 0.2484;

  a[0] = 1.4335e-4;
  a[1] = 1.2456e-4;
  a[2] = 8.1884e-5;
  a[3] = 5.9578e-5;
  a[4] = 4.9952e-5;
  a[5] = 4.2588e-5;

  b[0] = 3.3954e-9; 
  b[1] = 1.4612e-8;
  b[2] = 1.9394e-9;
  b[3] = 1.6707e-9;
  b[4] = 1.3893e-9;
  b[5] = 1.1963e-9;

  c[0] = -7.6364e-12;
  c[1] = -7.6345e-12;
  c[2] = -4.3618e-12;
  c[3] = -3.0058e-12;
  c[4] = -2.5331e-12;
  c[5] = -2.1709e-12;

  d[0] = -1.1950e-15;
  d[1] = -1.7511E-15;
  d[2] = -6.8260e-16;
  d[3] = -5.1925e-16;
  d[4] = -4.3896e-16;
  d[5] = -3.8282e-16;
 
  tstar = teff - 5700;

  for (i=0;i<6;i++) {
    seff[i] = seffsun[i] + a[i]*tstar + b[i]*tstar*tstar + 
c[i]*pow(tstar,3) + d[i]*pow(tstar,4);
   lim[i] = pow(l/seff[i],0.5);
  }
}

int main(int argc,char *argv[]) {
  int i;
  double dLum,dTeff,lim[6];

  if (argc != 5) {
    fprintf(stderr,"ERROR: Incorrect number of arguments. Usage: %s -l luminosity[solar] -t Teff[K].\n",argv[0]);
    exit(1);
  }

  dLum=0;
  dTeff=0;

  /* Check for flags */
  /* This really needs to be treated better */
  for (i=1;i<argc;i++) {
    if (memcmp(argv[i],"-l",2) == 0) 
      dLum = atof(argv[i+1]);
    if (memcmp(argv[i],"-t",2) == 0) 
      dTeff = atof(argv[i+1]);
    }

   if (dTeff == 0 || dLum == 0) {
     fprintf(stderr,"ERROR: Usage: %s -l luminosity[solar] -t Teff[K].\n",argv[0]);
     exit(1);
    }

  HabitableZone(dLum,dTeff,lim);

  printf("Recent Venus: %.4lf AU\n",lim[0]);
  printf("Runaway Greenhouse: %.4lf AU\n",lim[1]);
  printf("Moist Greenhouse: %.4lf AU\n",lim[2]);
  printf("Maximum Greenhouse: %.4lf AU\n",lim[3]);
  printf("Early Mars: %.4lf AU\n",lim[4]);

  printf("HZ: %.4lf -- %.4lf\n",(lim[0]+lim[2])/2,(lim[3]+lim[4])/2);
}