/* 
   INTENSITY_MEAN.cxx

   Version 1.0

   Calculate the intensity mean and the phase-weighted intensity mean. 
   Also calculate the difference between these two means.

   Reference: a) Saha & Hoessel, 1990, AJ 99:97
              b) Ferrarese et al., 1996, ApJ 464:568 

   Equations/algorithms:
   A) Intensity mean:
         1. m_av = -2.5 log10[ (1/N) x sum 10^(-0.4 m_i) ], i=1,..,N

   B) Phase-weighted intensity mean:
         1. Sort the phases.
         2. m_ph = -2.5 log10[ sum 0.5(phase_{i+1}-phase_{i-1})
	                x 10^(-0.4 m_i) ], i=1,..,N    
         3. Set phase_{0} = phase_{N} - 1
                phase_{N+1} = phase_{1} + 1

   Required sub-directories:
         1. /NGC_#, where # is the NGC number
	 2. /NGC_#/Vphase/
   
   Reqiured input files:
         1.  /NGC_#/period.#
	 2.  /NGC_#/Vphase/datav_p.#.i, where i is the i^th Cepheid

  Has not included the calculation/estimation of errors in the means. 

  Compile: g++ -o intmean intensity_mean.cxx

  Writte by C. Ngeow
  06/05/2003
*/

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

#define N 350

/* Sort the phase points */
void sort(int n, float phase[], float mag[], float sig[]);

void sort(int n, float phase[], float mag[], float sig[])
{
  int i, j;
  float a, b, c;

  for(j=2; j<=n; j++) {
    a = phase[j];
    b = mag[j];
    c = sig[j];

    i = j-1;
    
    while(i>0 && phase[i]>a)
      {
	phase[i+1] = phase[i];
	mag[i+1] = mag[i];
	sig[i+1] = sig[i];
	i--;
      }
    phase[i+1] = a;
    mag[i+1] = b;
    sig[i+1] = c;
  }
}

main()
{
  int ngc;

  cout<<" ** Intensity Mean & Phase-Weighted Intensity Mean ** \n\n";
  cout<<"  Input NGC : ";
  cin>>ngc;

  int Nceph;

  /* Get Nceph */

  char fper[30];

  sprintf(fper, "NGC_%d/period.%d", ngc, ngc);
  ifstream per(fper);

  int count = 0;
  float period;
  
  while(per>>period)
    count++;

  Nceph = count;

  /* ********************** */
  
  int n;
  float tp, tm, ts;
  float phase[N];
  float mag[N];
  float sig[N];

  float sum_mav, sum_mph;
  float m_av, m_ph;
  float delta;

  char fin[60];

  cout<<"\nV band results for NGC"<<ngc<<":\n";
  cout<<"Ceph\tNepoch\tm_av\tm_ph\tdelta\n";
  cout<<"------------------------------------------\n";

  for(int j=1; j<=Nceph; j++) {
    
    /* Input data and sort the phase values */

    sprintf(fin, "NGC_%d/Vphase/datav_p.%d.%d", ngc, ngc, j);
    ifstream in(fin);

    n = 0;
    while(in>>tp>>tm>>ts)
      {
	if(tm < 95.0)
	  {
	    phase[n+1] = tp;
	    mag[n+1] = tm;
	    sig[n+1] = ts;
	    n++;
	  }
	else
	  n = n;
      }
    
    /* Sort the data */
    
    sort(n, phase, mag, sig);

    /* Calculate intensity mean */

    sum_mav = 0.0;
    for(int i=1; i<=n; i++) {
      sum_mav += (1.0/float(n)) * pow(10.0, -0.4*mag[i]);
    }

    m_av = -2.5 * log10(sum_mav);

    /* Calculate phase-weighted intensity mean */

    sum_mph = 0.0;
   
    phase[0] = phase[n] - 1.0;
    phase[n+1] = phase[1] + 1.0;

    for(int i=1; i<=n; i++) 
      sum_mph += 0.5*(phase[i+1]-phase[i-1]) * pow(10.0, -0.4*mag[i]);
    
    m_ph = -2.5 * log10(sum_mph);

    /* Calculate delta */
    
    delta = fabs(m_av - m_ph);

    /* Output data */

    printf("%d\t%d\t%2.3f\t%2.3f\t%2.3f\n", j, n, m_av, m_ph, delta);

    in.close();
  }
}

#undef N