#include "star.hxx"

star::star(){  // constructor (set invalid values)
  icatnum = -999;
  ra_h = -999.;
  dec_deg = -999.;
  umag = -999.;
  bmag = -999.;
  vmag = -999.;
  rmag = -999.;
  u = -999.;
  v = -999.;
  ra_rad = -999.;
  dec_rad = -999.;
}

int star::readstar(FILE *fp, int verbose){ // read one line of the SKY2000 V2.0
                              // catalog and extract the interesting
                              // data
  int ira_hours, ira_min;
  int idec_degrees, idec_arcmin;
  float ra_sec, dec_arcsec;
  char catline[SKY2000LINELENGTH + 1];
  char *pos;
  char c2[3];
  char c3[4];
  char c6[7];
  char c7[8];
  char c8[9];

  strncpy(c2+2, "\0", 1);
  strncpy(c3+3, "\0", 1);
  strncpy(c6+6, "\0", 1);
  strncpy(c7+7, "\0", 1);
  strncpy(c8+8, "\0", 1);

  pos = catline;

  if( fgets( pos , SKY2000LINELENGTH + 1, fp) == NULL ){
    return(FALSE);
  }

  if(verbose > 2) fprintf(stdout, "%s\n", catline); 

  pos = catline + 27;


  strncpy(c8, pos, 8);
  sscanf(c8, "%d", &icatnum);

  pos = catline + 118;
  strncpy(c2, pos, 2);
  sscanf(c2, "%d", &ira_hours);

  pos = catline + 120;
  strncpy(c2, pos, 2);
  sscanf(c2, "%d", &ira_min);

  pos = catline + 122;
  strncpy(c7, pos, 7);
  sscanf(c7, "%f", &ra_sec);

  pos = catline + 129;
  strncpy(c3, pos, 3);
  if( c3[1] == ' ' ){
    c3[1] = '0';
  }
  if( c3[2] == ' ' ){
    c3[2] = '0';
  }
  sscanf(c3, "%d", &idec_degrees);

  pos = catline + 132;
  strncpy(c2, pos, 2);
  sscanf(c2, "%d", &idec_arcmin);

  pos = catline + 134;
  strncpy(c6, pos, 6);
  sscanf(c6, "%f", &dec_arcsec);

  pos = catline + 231;
  strncpy(c6, pos, 6);
  if (sscanf(c6, "%f", &vmag) < 1)
    vmag = -999.;;

  pos = catline + 251;
  strncpy(c6, pos, 6);
  if (sscanf(c6, "%f", &bmag) < 1)
    bmag = -999.;

  pos = catline + 271;
  strncpy(c6, pos, 6);
  if (sscanf(c6, "%f", &umag) < 1)
    umag = -999.;

  ra_h = ira_hours + ira_min/60. + ra_sec/3600.;
  dec_deg = idec_degrees + idec_arcmin/60. + dec_arcsec/3600.;
  ra_rad = ra_h * PI / 12.;
  dec_rad = dec_deg * PI /180.;

  if (verbose > 2) fprintf(stdout, "extracted: %d %d %d %f %d %d %f %f %f %f\n", 
		       icatnum, ira_hours, ira_min, ra_sec, 
		       idec_degrees, idec_arcmin, dec_arcsec, 
		       umag, bmag, vmag);
  
  return(TRUE);
}

int star::printstar(){ // write one star's parameters
  fprintf(stdout, "%d %f %f %f %f %f %f\n", 
 	  icatnum, ra_h, dec_deg, umag, bmag, vmag, rmag);
  return(0);
}

//----------------------------------------------------------------------------
// @name calcmissingmags
//
// @desc calculate the magnitudes for those wavebands in which no data
// @desc is available assuming a black body and using the V and B mags
//
//----------------------------------------------------------------------------

float star::calcmissingmags(int verbose) { // returns effective temperature; -1. = not possible 

  float temp;
  float tprime;
  float nu1_Hz, nu2_Hz, bflux, vflux, rflux, xmag;



  if(vmag > -100.){ // valid vmag  available
    if(bmag < -100.){ // no valid bmag  available
      cout << "Warning: star no. " <<  icatnum << 
	" has no Bmag measurement. Using Bmag = Vmag = "<<vmag<<"\n";
      bmag = vmag;
    }
  }
  else{
    cout << "Warning: star no. " <<  icatnum << 
	" has no Vmag measurement.\n";
    return(-1.);
  }
 
  // calculate the star temperature using approximation from
  // Kitchin, C.R., Astrophysical Techniques, 2nd ed., equ. 3.1.24

  if((bmag-vmag) > -0.2){
    temp = 8540. / ( (bmag-vmag) + 0.865 );
    if (verbose > 1) cout << "Star temperature from B-V: T = " << temp << "K\n";
  }
  else{
    temp = 12000.;
    if (verbose > 1) cout << "Star temperature from B-V: T > " << temp << "K\n";
  }

  // calculate an effective temperature for the Rmag calculation 
  // tprime = T * k / h

  nu1_Hz = LIGHTSPEED_mps/((VLMIN_nm+VLMAX_nm)/2.*1e-9);
  nu2_Hz = LIGHTSPEED_mps/((BLMAX_nm+BLMAX_nm)/2.*1e-9);
  vflux = pow(10.,-0.4*vmag-22.42);
  bflux = pow(10.,-0.4*bmag-22.42);
  tprime = (nu2_Hz - nu1_Hz) / ( log(vflux/bflux) - 3. * log(nu1_Hz/nu2_Hz) );
  if (verbose > 1) cout << "Blackbody T = " << tprime/1.38e-23*6.62e-34 << "\n";
     
     
  if( umag < -100. ){ // umag could not be read
    if (verbose) cout << "Warning: star no. " <<  icatnum << 
      " has no Umag measurement. Calculating it from its Vmag = " << vmag << "\n";
    if (verbose) cout << "         and Bmag = " << bmag << " assuming standard colour-colour-plot  ... ";

      if((bmag-vmag) > 1.4){
	umag =  bmag * 0.9;
      }
      else{
	if((bmag-vmag) > 0.5){
	  umag = -0.5 + 1.37 * (bmag - vmag) + bmag; 
	}
	else{
	  if((bmag-vmag) <= 0.){
	    umag = 4.07 * (bmag - vmag) + bmag;
	  }
	  else{
	    umag = 0.175 * (bmag - vmag) + bmag;
	  }
	}
      }

    if (verbose) cout << " result Umag = " << umag << "\n";
    
    if( umag < 5.0 ){
      cout << "Warning: star no. " <<  icatnum << " is bright (Vmag =" << vmag << ", Bmag = "
	   << bmag << ")\n         and has no Umag measurement. Estimated Umag is "<< umag <<"\n";
    }

  }
  else{ // umag available
    if (verbose > 1) {
      cout << "Test: star no. " <<  icatnum << 
	" has Umag = " << umag <<". Calculating it from its Vmag = " << vmag << "\n";
      cout << "         and Bmag = " << bmag << " assuming standard colour-colour-plot ...\n ";
      
      if((bmag-vmag) > 1.4){
	xmag =  bmag * 0.9;
      }
      else{
	if((bmag-vmag) > 0.5){
	  xmag = -0.5 + 1.37 * (bmag - vmag) + bmag; 
	}
	else{
	  if((bmag-vmag) <= 0.){
	    xmag = 4.07 * (bmag - vmag) + bmag;
	  }
	  else{
	    xmag = 0.175 * (bmag - vmag) + bmag;
	  }
	}
      }
      if (verbose > 2) 
	cout << "TEST " << umag <<" "<< xmag << " " << temp << " " << bmag << " " << vmag << "\n";

      cout << "          result Umag = " << xmag << "\n";
    }

  }

  if( rmag < -100. ){ // rmag not present (it's not part of the catalog)

    temp = temp * 1.38e-23 / 6.62e-34;     // * k / h

    rflux = vflux * pow((VLMIN_nm + VLMAX_nm)/(RLMIN_nm + RLMAX_nm),3.) *
      (exp(nu1_Hz/temp) - 1.) /
      (exp(LIGHTSPEED_mps/((RLMIN_nm+RLMAX_nm)/2.*1e-9)/temp) - 1.);

    rmag = (log10(rflux) + 22.42)/(-0.4);

  }

  return(tprime);
}
  

//----------------------------------------------------------------------------
// @name mag_nphot
//
// @desc translates magnitudes in number of photons, using log(flux)= -0.4*m+22.42
//
//----------------------------------------------------------------------------

int star::mag_nphot(int np[4], float inttime_s, float radius_m, int verbose) {
  
  float bflux, vflux, uflux, rflux; 
  float bintensity, vintensity, uintensity, rintensity;
  float unu1_Hz, unu2_Hz, bnu1_Hz, bnu2_Hz, vnu1_Hz, vnu2_Hz, rnu1_Hz, rnu2_Hz;

  //	The flux is given in Watt/m2*Hz.
  
  uflux = pow(10.,-0.4*umag-22.42);
  bflux = pow(10.,-0.4*bmag-22.42);
  vflux = pow(10.,-0.4*vmag-22.42);
  rflux = pow(10.,-0.4*rmag-22.42);

  if (verbose) cout<<"MAGS "<<umag<<" "<<bmag<<" "<<vmag<<" "<<rmag<<endl;
  
  unu1_Hz = LIGHTSPEED_mps/(ULMIN_nm*1e-9);
  unu2_Hz = LIGHTSPEED_mps/(ULMAX_nm*1e-9);
  bnu1_Hz = LIGHTSPEED_mps/(BLMIN_nm*1e-9);
  bnu2_Hz = LIGHTSPEED_mps/(BLMAX_nm*1e-9);
  vnu1_Hz = LIGHTSPEED_mps/(VLMIN_nm*1e-9);
  vnu2_Hz = LIGHTSPEED_mps/(VLMAX_nm*1e-9);
  rnu1_Hz = LIGHTSPEED_mps/(RLMIN_nm*1e-9);
  rnu2_Hz = LIGHTSPEED_mps/(RLMAX_nm*1e-9);

  // The intensity is given in number_of_photons/sec*m2. We obtain this units
  // because we multiply by the conversion factor 1Joule/s=h*c/((lambda1+lambda2)/2)
  
  
  uintensity = uflux*(unu1_Hz-unu2_Hz) * (ULMIN_nm+ULMAX_nm)*1e-9/2. /(PLANCK_si*LIGHTSPEED_mps);
  bintensity = bflux*(bnu1_Hz-bnu2_Hz) * (BLMIN_nm+BLMAX_nm)*1e-9/2. /(PLANCK_si*LIGHTSPEED_mps);
  vintensity = vflux*(vnu1_Hz-vnu2_Hz) * (VLMIN_nm+VLMAX_nm)*1e-9/2. /(PLANCK_si*LIGHTSPEED_mps);
  rintensity = rflux*(rnu1_Hz-rnu2_Hz) * (RLMIN_nm+RLMAX_nm)*1e-9/2. /(PLANCK_si*LIGHTSPEED_mps);

  np[0] = (int)(uintensity * radius_m * radius_m * PI * inttime_s);  
  np[1] = (int)(bintensity * radius_m * radius_m * PI * inttime_s);  
  np[2] = (int)(vintensity * radius_m * radius_m * PI * inttime_s);  
  np[3] = (int)(rintensity * radius_m * radius_m * PI * inttime_s);  

  if (verbose) cout<<"NPH "<<np[0]<<" "<<np[1]<<" "<<np[2]<<" " <<np[3]<<endl;

  numphot = np[0] + np[1] + np[2] + np[3];

  return(numphot);
  
}
