#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <libgen.h> // for basename()

int main(int argc,char **argv){
 if( argc<5 ){
  fprintf(stderr,"Usage: %s OBSERVED_MAG ERR BAND E(B-V)\nEample: %s 18.61 0.15 UVW1 0.011",argv[0],argv[0]);
  return 1;
 }
 double mag = atof(argv[1]);
 double mag_err = atof(argv[2]);
 double EBmV=atof(argv[4]);
 double Rv=3.1;
 double A,a,b;
 double c_m_sec = 299792458; /* speed of light m/sec */
 //
 double lambda_eff,nu,flux_convertion_factor,Zpt,flux,flux_max,flux_min; 
 //
 int band_id=0; // 0 - unknown band; 1 - known band
 /* 
      Set filter parameters in the fllowing lines.
 
            Zero_mag_flux * W/m^2/Hz->erg/sec/cm^2/Hz  */

 if( strlen(argv[3])<3 ){
  if( band_id==0 && 0==strncmp("V",argv[3], strlen(argv[3]) ) ){
   lambda_eff = 0.5402 * 1e-6; /* micron to meter */
   nu = c_m_sec/lambda_eff;
   flux_convertion_factor = 2.614e-16 * c_m_sec*1e10/(nu*nu);
   Zpt = 17.89;
   a = 1.0015; b = 0.0126;
   band_id=1;
  }
  if( band_id==0 && 0==strncmp("B",argv[3], strlen(argv[3]) ) ){
   lambda_eff = 0.4329 * 1e-6; /* micron to meter */
   nu = c_m_sec/lambda_eff;
   flux_convertion_factor = 1.472e-16 * c_m_sec*1e10/(nu*nu);
   Zpt = 19.11;
   a = 0.9994; b = 1.0171;
   band_id=1;
  } 
  if( band_id==0 && 0==strncmp("U",argv[3], strlen(argv[3]) ) ){
   lambda_eff = 0.3501 * 1e-6; /* micron to meter */
   nu = c_m_sec/lambda_eff;
   flux_convertion_factor = 1.63e-16 * c_m_sec*1e10/(nu*nu);
   Zpt = 18.34;
   a = 0.9226; b = 2.1019;
   band_id=1;
  }
 }
 if( band_id==0 && 0==strncmp("UVW1",argv[3], strlen(argv[3]) ) ){
  lambda_eff = 0.2634 * 1e-6; /* micron to meter */
  nu = c_m_sec/lambda_eff;
  //flux_convertion_factor = 4.3e-16 * c_m_sec*1e10/(nu*nu);
  flux_convertion_factor = 4.0e-16 * c_m_sec*1e10/(nu*nu);
  Zpt = 17.49;
  a = 0.4346; b = 5.3286;
  band_id=1;
 }
 if( band_id==0 && 0==strncmp("UVM2",argv[3], strlen(argv[3]) ) ){
  lambda_eff = 0.2231 * 1e-6; /* micron to meter */
  nu = c_m_sec/lambda_eff;
  //flux_convertion_factor = 7.5e-16 * c_m_sec*1e10/(nu*nu);
  flux_convertion_factor = 8.5e-16 * c_m_sec*1e10/(nu*nu);
  Zpt = 16.82;
  a = 0.0773; b = 9.1784;
  band_id=1;
 }
 if( band_id==0 && 0==strncmp("UVW2",argv[3], strlen(argv[3]) ) ){
  lambda_eff = 0.2030 * 1e-6; /* micron to meter */
  nu = c_m_sec/lambda_eff;
  //flux_convertion_factor = 6.0e-16 * c_m_sec*1e10/(nu*nu);
  flux_convertion_factor = 6.2e-16 * c_m_sec*1e10/(nu*nu);
  Zpt = 17.35;
  a = -0.0581; b = 8.4402;
  band_id=1;
 }
 
 if( band_id==0 ){
  fprintf(stderr,"ERROR: unknown band %s\n",argv[3]);
  return 1;
 }
 
 A=EBmV*(a*Rv+b); // extinction in this band
 mag=mag-A; // apply extinction correction
 
 // --------------------------------------------------------- //
 
 if( 0==strcmp("UVOT2ergcm2sA_wavelength_Poole2008",basename(argv[0])) ){
  // https://en.wikipedia.org/wiki/Planck%27s_law
  flux_convertion_factor=flux_convertion_factor*(nu*nu)/(c_m_sec*1e10);
 }
 
 // pow(10,0.4*(Zpt-mag)) is the XRT count rate
 // We need to multiply it by the flux conversion factor to get the flux
 flux = flux_convertion_factor * pow(10,0.4*(Zpt-mag));
 flux_max = flux_convertion_factor * pow(10,0.4*(Zpt-mag-mag_err));
 flux_min = flux_convertion_factor * pow(10,0.4*(Zpt-mag+mag_err));

 if( 0==strcmp("UVOT2ergcm2sA_wavelength_Poole2008",basename(argv[0])) ){
  fprintf(stdout,"%7.2lf  %lg  %lg\n", 1e10*c_m_sec/nu, flux, fabs(flux_max-flux_min)/2.0  );
 }
 else{
  fprintf(stdout,"%lg  %lg  %lg\n",nu,flux*pow(10,23)*pow(10,6),( fabs(flux_max-flux_min)/2.0 )*pow(10,23)*pow(10,6));
 }
  
 return 0;
}