/* CYCLOCODE: Period Searching & Lightcurve Fitting Program */
/* Ver. 3.6; B. Dermawan */
/*
This program consists of two parts, which can be employed separately or simultaneously, 
depends on user's choice. User inputs some parameters when employing this program.
Part 1 evaluates a period of the lightcurve data using the PDM and/or the Spectral Analysis methods.
   The period is calculated in unit day (frequency of 1/day).
Part 2 develops appropriate fitting curve(s) of the lightcurve. Base-level of a fitting curve can be
   both at about the mean level of the magnitude data and at zero magnitude by shifting of delta magnitude.
   The fitting curve/data for the latter is called "normalized" magnitude (normag).
The input file (*.???) must have 3 columns and time-increasing sorted: 
    Column 1 -> Time in JD/MJD (or its fraction) relative to epoch,
    Column 2 -> Magnitude (or relative magnitude),
    Column 3 -> Error of magnitude (or relative magnitude).
The output files: Part 1 -> *.pdm: PDM-based period result
                         -> *.spc: Spectral-based period result
                  Part 2 -> *.pha: Rotational-phased data
                         -> *.roc: Rotational-phased fitting curve
                         -> *.fuc: Time-span fitting (continuous) curve 
                         -> *.fud: Time-span fitting data points at the same 'Time' data
This code is open, freely improved and distributed.
*/
#include <stdio.h> 
#include <stdlib.h> 
#include <math.h> 
#include <string.h>
#include <ctype.h>
#define NR_END 1                                /*Constants*/
#define FREE_ARG char*                          
#define TWOPID 6.2831853071795865  	            
#define TINY 1.0e-20                            
#define DATP 80                                
#define OFF 0.25                                
#define SWAP(a,b) {swap=(a);(a)=(b);(b)=swap;}  /*Defined subroutine*/

int *ivector(long nl, long nh);                                                       /*Subroutine Names*/
float *vector(long nl, long nh);                                                                  
float **matrix(long nrl, long nrh, long ncl, long nch);                                            
void free_vector(float *v, long nl, long nh);                                                      
void free_ivector(int *v, long nl, long nh);                                                       
void free_matrix(float **m, long nrl, long nrh, long ncl, long nch);                               

float SQR(float);                                                                                  
void PDM(float time[], float mag[], FILE *pdmPtr, float *pdmper, int n);                           
void SpAn(float time[], float mag[], FILE *spcPtr, float *spanper, int n);                         
void Sort(float *vart, float *varc1, float *varc2, float *varc3, int ndat);                        
void LinFit(float x[], float y[], float sig[], float x0, int ndat, float *b, int ma, float *afunc, 
            float *chisq, float Period, void (*)(float, float [], float, int, float));             
void FouHarm(float x, float *afunc, float x0, int ma, float Period);                               
void FitCurve(float x[], float *y, int datarange, float period, float x0, float *b, int ma,        
              float *afunc, void (*FouHarm)(float, float [], float, int, float));                  

main()             /* ### MAIN PROGRAM ### */
{ 
   char c,choice[3],fin[30],fpdm[30],fspc[30],fpha[30],froc[30],ffuc[30],ffud[30];              /*Variables*/
   char pdmf[]="pdm",spcf[]="spc",phaf[]="pha",rocf[]="roc",fucf[]="fuc",fudf[]="fud";             
   int ce=3,cm,cp,ndat=1,f,ord,lo,up,datr,ma,i,j,k;                                         
   float *ph,*mag,*JD,*sig,*b,*afunc,*xd,*yd,*chisq,**phamatr,**rocmatr,**fucmatr,**fudmatr;       
   float aa,bb,cc,pdmper,spanper,fullper,Period,x0,stap,endp;                                     
   FILE *datPtr,*pdmPtr,*spcPtr,*phaPtr,*rocPtr,*fucPtr,*fudPtr;                                   

   printf("\nThis program evaluates: 1. Periodogram(s); 2. Synthetic Lightcurve(s); 3. Both\n");
   printf("Please input your choice: "); 
   j=0; while((c=getchar()) != '\n') choice[j++]=c; choice[j]='\0'; sscanf(choice,"%d",&ce);
   if ((ce != 1) && (ce != 2) && (ce != 3)) ce = 3;
   
   j = 0; printf("\nInput File (*.\?\?\?): ");                         /*Filename Assemblages*/
   while((c=getchar()) != '\n') {                                      
   	  fin[j]=fpdm[j]=fspc[j]=fpha[j]=froc[j]=ffuc[j]=ffud[j]=c; j++;   
   }                                                                   
   fin[j]=fpdm[j]=fspc[j]=fpha[j]=froc[j]=ffuc[j]=ffud[j]='\0';        
   for(k=3;k>=1;k--) {                                                 
      fpdm[j-k]=pdmf[3-k]; fspc[j-k]=spcf[3-k];                        
      fpha[j-k]=phaf[3-k]; froc[j-k]=rocf[3-k];                        
      ffuc[j-k]=fucf[3-k]; ffud[j-k]=fudf[3-k];                       
   }                                                                   

   if ((datPtr = fopen(fin,"r")) == NULL) {             
      printf("File could not be opened\n"); return 0;	
   }                                                    
   else {                                               
      fscanf(datPtr,"%f %f %f", &aa, &bb, &cc);                        /*Checking number of data*/
      while (!feof(datPtr)) {
         fscanf(datPtr,"%f %f %f", &aa, &bb, &cc); ndat++;   	
      }
   }
   printf("Number of data = %3d\n",ndat-=1);
   
   JD=vector(1,ndat); mag=vector(1,ndat); sig=vector(1,ndat); datPtr = fopen(fin,"r");
   for (i=1;i<=ndat;i++) fscanf(datPtr, "%f %f %f", &JD[i], &mag[i], &sig[i]);             /*Reading the input file*/

   if (ce != 1) {
      if (JD[ndat]-JD[1] > 0.6) {                          /*Checking the time-pan, single night or not*/
         printf ("Base-level of the mag. data is zero mag. [y/Y (default) or n/N]: ");
         j=0; while((c=getchar()) != '\n') choice[j++]=c; choice[j]='\0'; 
         if ((choice[0] == 'n') || (choice[0] == 'N')) f=2; else f=1;
      } else f=2;
   }

   if ((ce == 1) || (ce == 3)) { 
      printf("\nPERIODOGRAM(S)\n");
      printf("==============\n");
      printf("   Method: [1] PDM; [2] Spectral Analysis (default); [3] Both: ");  
      j=0; while((c=getchar()) != '\n') choice[j++]=c; choice[j]='\0'; sscanf(choice,"%d",&cm); printf("\n");
      switch (cm) {
         case 1: pdmPtr=fopen(fpdm,"w"); PDM(JD,mag,pdmPtr,&pdmper,ndat); fullper=pdmper; 
                 if (ce == 1) {
                    fclose(pdmPtr); fclose(datPtr); free_vector(JD,1,ndat); free_vector(mag,1,ndat); free_vector(sig,1,ndat); 
                    printf("\n"); return 0;
                 }
                 break;
         case 2: default: spcPtr=fopen(fspc,"w"); SpAn(JD,mag,spcPtr,&spanper,ndat); fullper=spanper; 
                 if (ce == 1) {
                    fclose(spcPtr); free_vector(JD,1,ndat); free_vector(mag,1,ndat); free_vector(sig,1,ndat); fclose(datPtr); 
                    printf("\n"); return 0;
                 }
                 break;
         case 3: pdmPtr=fopen(fpdm,"w"); PDM(JD,mag,pdmPtr,&pdmper,ndat); 
                 spcPtr=fopen(fspc,"w"); SpAn(JD,mag,spcPtr,&spanper,ndat);
                 if (ce == 1) {
                    fclose(pdmPtr); fclose(spcPtr); free_vector(JD,1,ndat); free_vector(mag,1,ndat); free_vector(sig,1,ndat); 
                    fclose(datPtr); printf("\n"); return 0;
                 }
                 printf("\nThe obtained period used for developing synthetic lightcurve(s)\n");
                 printf("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");
                 printf("     [1] PDM; [2] Spectral Analysis (default): "); 
                 i=0; while((c=getchar()) != '\n') choice[i++]=c; choice[i]='\0'; sscanf(choice,"%d",&cp);
                 if (cp == 1) fullper=pdmper; else fullper=spanper;
                 break;
      }
   } 

   printf("\nSYNTHETIC LIGHTCURVE(S)\n");
   printf("=======================\n");
   if (ce == 2) { printf ("   Known period (day): "); scanf("%f",&fullper); printf("\n"); }
   printf("   Least-squares Fitting:\n");
   printf("     Range of Fourier-harmonics' order: Lower = "); scanf("%d", &lo);
   printf("                                        Upper = "); scanf("%d", &up);
   printf("     Order   Base   Chi-Square\n");
   printf("     -------------------------\n");


   ph=vector(1,ndat); chisq=vector(1,up-lo+1); 
   phaPtr=fopen(fpha,"w"); rocPtr=fopen(froc,"w"); fucPtr=fopen(ffuc,"w"); fudPtr=fopen(ffud,"w");
   fprintf(rocPtr,"            "); fprintf(phaPtr,"                                    ");
   fprintf(fucPtr,"                     "); fprintf(fudPtr,"                     ");
   
   rocmatr=matrix(1,1+DATP,1,2*(up-lo+1)+1); phamatr=matrix(1,ndat,1,up-lo+5);
   fucmatr=matrix(1,1+ceil(DATP*(JD[ndat]-JD[1]+2.0*OFF*fullper)/fullper),1,2*(up-lo+2)); 
   fudmatr=matrix(1,ndat,1,2*(up-lo+2)); 
   for (i=1;i<=ndat;i++) ph[i]=fmod(JD[i],fullper)/fullper;            /*Rotational-phased data*/		

   for (ord=lo;ord<=up;ord++) { 
      ma=ord*2+1; afunc=vector(1,ma); b=vector(1,ma);                  /*Basis function allocation*/
      
      Sort(ph,JD,mag,sig,ndat); Period=1.0; x0=(ph[1]+ph[ndat])/2.0;     
      LinFit(ph,mag,sig,x0,ndat,b,ma,afunc,chisq,Period,FouHarm);	   /*Rotational-phased fitting task*/
      printf("%8d %9.4f %11.3E\n",ord,b[ma],*chisq);

      for (i=1;i<=4;i++) {              /*Storing the results into the output matrices*/
         switch (i) {
            case 1: stap=0.0; endp=1.0; datr=1+(int)(DATP*(endp-stap)/Period); 
                    xd=vector(1,datr); yd=vector(1,datr); for (k=1;k<=datr;k++) xd[k]=stap+(k-1)*Period/DATP;
                    if (ord == lo) {
                       for (k=1;k<=datr;k++) rocmatr[k][1]=stap+(k-1)*Period/DATP;
                       for (k=1;k<=ndat;k++) {
                          phamatr[k][1]=ph[k]; phamatr[k][2]=JD[k]; phamatr[k][3]=mag[k]; phamatr[k][4]=sig[k];
                       }
                    }
                    break;
            case 2: datr=ndat; xd=vector(1,datr); yd=vector(1,datr); for (k=1;k<=datr;k++) xd[k]=ph[k]; break;
            case 3: Sort(JD,ph,mag,sig,ndat); Period=fullper; x0=Period/2.0; 
                    stap=JD[1]-OFF*Period; endp=JD[ndat]+OFF*Period; 
                    datr=1+(int)(DATP*(endp-stap)/Period); xd=vector(1,datr); yd=vector(1,datr); 
                    for (k=1;k<=datr;k++) xd[k]=fmod(stap+(k-1)*Period/DATP,Period); 
                    if (ord == lo) {
                       for (k=1;k<=datr;k++) { fucmatr[k][1]=stap+(k-1)*Period/DATP; fucmatr[k][2]=xd[k]/Period; }
                       for (k=1;k<=ndat;k++) { fudmatr[k][1]=JD[k]; fudmatr[k][2]=ph[k]; }
                    } 
                    break;
            case 4: datr=ndat; xd=vector(1,datr); yd=vector(1,datr); for (k=1;k<=datr;k++) xd[k]=ph[k]*Period; break;
         }
         FitCurve(xd,yd,datr,Period,x0,b,ma,afunc,FouHarm); 
         for (k=1;k<=datr;k++) { 
            switch (i) {
               case 1: rocmatr[k][2*(ord-lo+1)]=yd[k]; rocmatr[k][2*(ord-lo+1)+1]=yd[k]-b[ma]; break;
               case 2: phamatr[k][ord-lo+5]=mag[k]-b[ma]; break;
               case 3: fucmatr[k][2*(ord-lo+2)-1]=yd[k]; fucmatr[k][2*(ord-lo+2)]=yd[k]-b[ma]; break;
               case 4: fudmatr[k][2*(ord-lo+2)-1]=yd[k]; fudmatr[k][2*(ord-lo+2)]=yd[k]-b[ma]; break;
            }
         }
         free_vector(xd,1,datr); free_vector(yd,1,datr); 
      }
      free_vector(b,1,ma); free_vector(afunc,1,ma); 
      fprintf(rocPtr," %2d    %8.4f   ",ord,b[ma]); fprintf(phaPtr,"  %8.4f",b[ma]); 
      fprintf(fucPtr," %2d    %8.4f   ",ord,b[ma]); fprintf(fudPtr," %2d    %8.4f   ",ord,b[ma]);
   }
                                        /*Constructing header of the output files*/
   fprintf(rocPtr,"\n"); fprintf(phaPtr,"\n"); fprintf(fucPtr,"\n"); fprintf(fudPtr,"\n");
   fprintf(rocPtr,"  Phase  "); fprintf(phaPtr,"  Phase      JD       mag      err  "); 
   fprintf(fucPtr,"    JD     Phase  "); fprintf(fudPtr,"    JD     Phase  ");
   for (ord=lo;ord<=up;ord++) {         /*Constructing header of the output files*/
      if (f == 1) {
         fprintf(rocPtr,"    mag      mag' "); fprintf(phaPtr,"   %2d%s",ord,":mag'");
         fprintf(fucPtr,"    mag      mag' "); fprintf(fudPtr,"    mag      mag' ");
      } else {
         fprintf(rocPtr,"    mag    normag "); fprintf(phaPtr," %2d%s",ord,":normag");
         fprintf(fucPtr,"    mag    normag "); fprintf(fudPtr,"    mag    normag ");
      }
   }
   fprintf(rocPtr,"\n"); fprintf(phaPtr,"\n"); fprintf(fucPtr,"\n"); fprintf(fudPtr,"\n");

   for (i=1;i<=ndat;i++) {              /*Arranging of the result matrices for output files*/
   	  for (j=1;j<=up-lo+5;j++) 
   	     if (j <= 2) fprintf(phaPtr,"%8.4f ",phamatr[i][j]);
   	     else if (j <= 4) fprintf(phaPtr,"%8.3f ",phamatr[i][j]); 
   	          else fprintf(phaPtr,"%9.3f ",phamatr[i][j]); 
   	  for (k=1;k<=2*(up-lo+2);k++) 
   	     if (k > 2) fprintf(fudPtr,"%8.3f ",fudmatr[i][k]); else fprintf(fudPtr,"%8.4f ",fudmatr[i][k]);
   	  fprintf(phaPtr,"\n"); fprintf(fudPtr,"\n");
   }
   for (i=1;i<=1+DATP;i++) {
   	  for (j=1;j<=2*(up-lo+1)+1;j++) 
   	     if (j > 1) fprintf(rocPtr,"%8.3f ",rocmatr[i][j]); else fprintf(rocPtr,"%8.4f ",rocmatr[i][j]);
   	  fprintf(rocPtr,"\n");
   }
   for (i=1;i<=1+(int)(DATP*(JD[ndat]-JD[1]+2.0*OFF*fullper)/fullper);i++) {
   	  for (k=1;k<=2*(up-lo+2);k++) 
   	     if (k > 2 ) fprintf(fucPtr,"%8.3f ",fucmatr[i][k]); else fprintf(fucPtr,"%8.4f ",fucmatr[i][k]);
   	  fprintf(fucPtr,"\n");
   }

   printf("\n"); free_matrix(rocmatr,1,1+DATP,1,f*(up-lo+1)+1); free_matrix(phamatr,1,ndat,1,up-lo+5); 
   free_matrix(fucmatr,1,1+ceil(DATP*(JD[ndat]-JD[1]+2.0*OFF*fullper)/fullper),1,2*(up-lo+2));
   free_matrix(fudmatr,1,ndat,1,2*(up-lo+2));
   fclose(datPtr); fclose(phaPtr); fclose(rocPtr); fclose(fucPtr); fclose(fudPtr);
   free_vector(ph,1,ndat); free_vector(mag,1,ndat); free_vector(JD,1,ndat); 
   free_vector(sig,1,ndat); free_vector(chisq,1,up-lo+1);
   
   return 0; 
} 

/* ### SUBROUTINES ### */

void Sort(float *varb, float *varm1, float *varm2, float *varm3, int n)
/* The last 3 columns (varm1, varm2, varm3) are sorted referring to the first column (varb) */
{ 
   int j,k;
   float swap; 
 
   for (k=1;k<=n-1;k++)  
      for (j=1;j<=n-1;j++) 
         if (varb[j] > varb[j+1]) { 
            SWAP(varb[j],varb[j+1]); SWAP(varm1[j],varm1[j+1]);
            SWAP(varm2[j],varm2[j+1]); SWAP(varm3[j],varm3[j+1]); 
         } 
}

void PDM(float time[], float magn[], FILE *pdmPtr, float *pdmper, int n)
/*Trial periods are examined from 0.001 d (~1.5 min) to 10 d. The minimum of Theta Statistics
means the most significance period. The routine returns a period (pdmper) that could be single-
or double-peak lightcurve.*/
{
   float Theta(float mg[], int n);
   
   int i,j,np;
   float *t,*m,*mm,*pdmphase,*Per,*Th,Pmin,Pmax,Pstep,Thmin=100.0,Prmin;
   
   printf("   - PDM:\n");
   t=vector(1,n); m=vector(1,n); mm=vector(1,n); pdmphase=vector(1,n);
   for (i=1;i<=n;i++) {
      t[i]=time[i]; m[i]=magn[i]; mm[i]=magn[i]+10.0;
   }
   Pmin=0.001; Pmax=10.0; Pstep=0.0003*(Pmin+Pmax)/2.0; np=1+ceil(Pmax-Pmin)/Pstep;
   Per=vector(1,np); Per[1]=Pmin; Th=vector(1,np);

   fprintf(pdmPtr,"Freq. (1/d)   Per. (d)   Theta\n");
   for (i=1;i<=np;i++) {
      for (j=1;j<=n;j++) pdmphase[j]=fmod(t[j],Per[i])/Per[i];
      Sort(pdmphase,t,mm,m,n);
      
      Th[i]=Theta(mm,n);
      if (Th[i] < Thmin) {
         Thmin=Th[i]; Prmin=Per[i];
      }
      Per[i+1]=Per[i]+Pstep;
   }
   fprintf(pdmPtr,"  Minimum:\n"); fprintf(pdmPtr,"%10.5f %10.4f %8.3f\n\n", 1.0/Prmin, Prmin, Thmin);
   for (i=1;i<=np;i++) {
      fprintf(pdmPtr,"%10.5f %10.4f %8.3f\n", 1.0/Per[i], Per[i], Th[i]);   /*Printing & Saving*/
   }
   
   *pdmper=Prmin; fclose(pdmPtr);
   printf("            Period = %8.4f %2s\n", Prmin, "day");
   printf("       Theta Stat. = %7.3f\n", Thmin);
   printf("-------------------------------\n");
   free_vector(Per,1,np); free_vector(Th,1,np);
   free_vector(t,1,n); free_vector(m,1,n); free_vector(mm,1,n); free_vector(pdmphase,1,n);
}

float Theta(float mg[], int n)
/*This Theta Statistics algorithm examines phase dispersion (Nsum/Dsum) of the data.*/ 
{
   int i;
   float Msum=0.0,Mmean=0.0,Nsum=0.0,Dsum=0.0;
 
   for (i=1;i<=n;i++) Msum+=mg[i];
   Mmean=Msum/n;
   for (i=1;i<=n-1;i++) {
      Nsum+=SQR(mg[i]-mg[i+1]); Dsum+=SQR(mg[i]-Mmean);
   }   
   return Nsum/Dsum;
}

void SpAn(float time[], float mag[], FILE *spcPtr, float *spanper, int n)
/*Oversampling parameter factor (ofac, >= 4) & high freq. (hifaq) are related as np = ofac*hifac*n/2,
where np is number of different freq. for n data points. User inputs the ofac and hifac values, and
select whether single- or double-peak lightcurve. The obtained period (spanper, longer than the Nyquist
interval) should have the highest spectral power with significance level of >= 99%. Otherwise, the data
may not have a periodic signal.*/
{
   void Spec(float x[], float y[], int n, float ofac, float hifac, float px[], float py[],
             int np, int *nout, int *jmax, float *prob, int *lim);
                
   char ch,choice[5];
   int sd=2,np,j,k,nout,jmax,*lim,peak;
   float ofac=4,hifac=15,*px,*py,prob,siglev;
   
   printf("   - Spectral Analysis:\n");
   printf("       Values of ofac & hifac ->  ofac = "); 
   j=0; while((ch=getchar()) != '\n') choice[j++]=ch; choice[j]='\0'; sscanf(choice,"%f",&ofac);	
   printf("                                 hifac = "); 
   j=0; while((ch=getchar()) != '\n') choice[j++]=ch; choice[j]='\0'; sscanf(choice,"%f",&hifac);     	

   printf("       Lightcurve peak: [Single => 1] or [Double => 2 (default)]: "); 
   j=0; while((ch=getchar()) != '\n') choice[j++]=ch; choice[j]='\0'; sscanf(choice,"%d",&sd); 
   if (sd == 1) peak=1; else peak=2;

   np=(int)(0.5*ofac*hifac*n)+1; px=vector(1,np); py=vector(1,np); lim=ivector(1,1);
   Spec(time,mag,n,ofac,hifac,px,py,np,&nout,&jmax,&prob,lim); 
   
   fprintf(spcPtr,"Freq. (1/d)   Per. (d)   Power  Sig. (%%)\n");         	/*Printing & Saving*/
   if (peak == 1) fprintf(spcPtr,"  Maximum:    (S-peak)\n");
   else fprintf(spcPtr,"  Maximum:    (D-peak)\n");
   fprintf(spcPtr,"%10.5f %10.4f %8.3f %8.2f\n\n",px[jmax]/peak,peak/px[jmax],py[jmax],(1.0-prob)*100.0);

   for (k=1;k<=*lim;k++) {
      siglev=1.0-hifac*n*exp(-py[k]);
      if (siglev < 0.99) siglev=pow(1.0-exp(-py[k]),hifac*n);
      fprintf(spcPtr,"%10.5f %10.4f %8.3f %8.2f\n", px[k]/peak, peak/px[k], py[k], siglev*100.0);
   }
   *spanper=peak/px[jmax];
   printf("            Period = %8.4f %2s\n", *spanper, "day");
   printf("        Sig. Level = %6.2f %3s\n", (1.0-prob)*100.0, "%");
   printf("-------------------------------\n");
   fclose(spcPtr); free_vector(px,1,np); free_vector(py,1,np); free_ivector(lim,1,1);
}

void FitCurve(float x[], float *y, int datarange, float period, float x0, float *b, int ma, 
              float *afunc, void (*FouHarm)(float, float [], float, int, float))
{
   int i,j;

   for (i=1;i<=datarange;i++) { 
      (*FouHarm)(x[i],afunc,x0,ma,period);      
      for (y[i]=0.0,j=1;j<=ma;j++) y[i] += b[j]*afunc[j];
   }
}

void Spec(float x[], float y[], int n, float ofac, float hifac, float px[], float py[],
          int np, int *nout, int *jmax, float *prob, int *lim) 
/*Given n data points with abscissas x[1..n] (which need not be equally spaced) and ordinates 
y[1..n], and given a desired oversampling factor ofac (a typical value being 4 or larger), 
this routine fills array px[1..np] with an increasing sequence of frequencies (not angular 
 frequencies) up to hifac times the "average" Nyquist frequency, and fills array py[1..np] 
with the values of the Lomb normalized periodogram at those frequencies. The arrays x and y 
are not altered. np, the dimension of px and py, must be large enough to contain the output, 
or an error results. The routine also returns jmax such that py[jmax] is the maximum element 
in py and prob an estimate of the significance of that maximum against the hypothesis of 
random noise. A small value of prob indicates that a significant periodic signal is present.*/ 
/*Small modifications from the original (Numerical Recipes in C)*/
{ 
   double *dvector(long nl, long nh);
   void free_dvector(double *v, long nl, long nh);
   void avevar(float data[], unsigned long n, float *ave, float *var);
   
   int i,j; 
   float ave,c,cc,cwtau,effm,expy,pnow,pymax,s,ss,sumc,sumcy,sums,sumsh, 
         sumsy,swtau,var,wtau,xave,xdif,dxmin=10.0,yy; 
   double arg,wtemp,*wi,*wpi,*wpr,*wr; 
    
   wi=dvector(1,n);  wpi=dvector(1,n);  wpr=dvector(1,n);  wr=dvector(1,n);
   
   *nout=(int)(0.5*ofac*hifac*n); 
   if (*nout > np) {	/*nrerror("Output arrays too short in period");*/
      printf("Output arrays too short in period\n"); exit(1); 
   }  
   avevar(y,n,&ave,&var);               /*Get mean and variance of the input data.*/ 
   if (var == 0.0) {	/*nrerror("Zero variance in period");*/
      printf("Zero variance in period\n"); exit(1);
   }  
   xdif=x[n]-x[1]; xave=0.5*(x[n]-x[1]);
   pymax=0.0;
   pnow=1.0/(xdif*ofac);                /*Starting frequency.*/ 
   for (j=1;j<=n;j++) {                 /*Initialize values for the trigonometric recurrences*/ 
      arg=TWOPID*((x[j]-xave)*pnow);    /*at each data point.The recurrences are done in double precision.*/ 
      wpr[j]=-2.0*SQR(sin(0.5*arg)); wpi[j]=sin(arg); wr[j]=cos(arg); wi[j]=wpi[j];

      if (j <= n-1) if (x[j+1]-x[j] < dxmin) dxmin=x[j+1]-x[j];    /*Find the minimum interval of the data*/
   } 
      
   for (i=1;i<=(*nout);i++) {           /*Main loop over the frequencies to be evaluated*/
      px[i]=pnow;                         
      sumsh=sumc=0.0;                   /*First,loop over the data to get tau and related*/ 
      for (j=1;j<=n;j++) {              /*quantities.*/ 
         c=wr[j]; s=wi[j]; 
         sumsh += s*c; sumc += (c-s)*(c+s); 
      } 
      wtau=0.5*atan2(2.0*sumsh,sumc); swtau=sin(wtau); cwtau=cos(wtau); 
      sums=sumc=sumsy=sumcy=0.0;       	/*Then,loop over the data again to get the*/ 
      for (j=1;j<=n;j++) {             	/*periodogram value.*/ 
         s=wi[j]; c=wr[j]; 
         ss=s*cwtau-c*swtau; cc=c*cwtau+s*swtau; 
         sums += ss*ss; sumc += cc*cc; 
         yy=y[j]-ave; 
         sumsy += yy*ss; sumcy += yy*cc; 
         wr[j]=((wtemp=wr[j])*wpr[j]-wi[j]*wpi[j])+wr[j];   /*Update the trigonometric recurrences.*/ 
         wi[j]=(wi[j]*wpr[j]+wtemp*wpi[j])+wi[j];  
      } 
      py[i]=0.5*(sumcy*sumcy/sumc+sumsy*sumsy/sums)/var; 
      if ((px[i] < 1.0/(2.0*dxmin)) && (py[i] >= pymax)) pymax=py[(*jmax=i)]; 
      *lim = i; 
      pnow += 1.0/(ofac*xdif);                             	/*The next frequency.*/ 
       
   } 
   expy=exp(-pymax); effm=2.0*(*nout)/ofac; 	/*Evaluate statistical significance of the maximum.*/ 
   *prob=effm*expy;  if (*prob > 0.01) *prob=1.0-pow(1.0-expy,effm); 

   free_dvector(wr,1,n); free_dvector(wpr,1,n); free_dvector(wpi,1,n); free_dvector(wi,1,n);
} 
 
void avevar(float data[], unsigned long n, float *ave, float *var) 
/*Given array data[1..n] returns its mean as ave and its variance as var.*/
/*Subroutine from Numerical Recipes in C*/ 
{ 
   unsigned long j; 
   float s,ep; 
 
   for (*ave=0.0,j=1;j<=n;j++) *ave += data[j]; 
   *ave /= n; *var=ep=0.0; 
   for (j=1;j<=n;j++) { 
      s=data[j]-(*ave); ep += s; *var += s*s; 
   } 
   *var=(*var-ep*ep/n)/(n-1);              	/*Corrected two-pass formula (14.1.8).*/ 
} 

void FouHarm(float x, float *afunc, float x0, int ma, float Period)      /* Fourier Harmonics function */
/*This routine returns the basis (afunc[k]) of Fourier Harmonics function for order ma*/
{ 
   int k;
    
   for (k=1;k<=ma-1;k++)            
      if (k <= (ma-1)/2) 
         afunc[k] = sin(k*TWOPID*(x-x0)/Period); 
      else afunc[k] = cos((k-(ma-1)/2)*TWOPID*(x-x0)/Period); 
   afunc[ma] = 1.0;          /*Phase angle effect (ma -> ma-1): afunc[ma-1] = TWOPID*(x-x0)/Period*/
}  

void LinFit(float x[], float y[], float sig[], float x0, int ndat, float *b, int ma, float *afunc,
            float *chisq, float Period, void (*FouHarm)(float, float [], float, int, float))
/*Given a set of data points x[1..ndat], y[1..ndat] with individual standard deviations 
sig[1..ndat], use Chi2 minimization to fit for some or all of the coefficients b[1..ma] of 
a function that depends linearly on b, y = Sigma_i b_i x afunc_i(x). The array ia[1..ma] indicates
by nonzero entries those components of b that should be fitted for, and by zero entries those
components that should be held fixed at their input values. The program returns values for b[1..ma],
Chi2 = chisq. (Parameters held fixed will return zero covariances.) The user supplies a routine
FouHarm(X,afunc,ma) that returns the ma basis functions evaluated at x = X in the array afunc[1..ma].*/ 
/* Small modifications from the original (Numerical Recipes in C) */
{ 
   void covsrt(float **covar, int ma, int ia[], int mfit); 
   void ludcmp(float **a, int n, int *indx, float *d); 
   void lubksb(float **a, int n, int *indx, float *b); 

   int *ia,*indx,ii,i,j,k,l,m,mfit=0; 
   float **covar,**beta,*a,*d,ym,wt,sum,sig2i; 
 
   covar=matrix(1,ma,1,ma); ia=ivector(1,ma); for (i=1;i<=ma;i++) ia[i]=1;
   beta=matrix(1,ma,1,1); a=vector(1,ma); d=vector(1,ma); indx=ivector(1,ma); 
   
   for (j=1;j<=ma;j++) 
      if (ia[j]) mfit++; 
   if (mfit == 0) {	/*nrerror ("LinFit: no parameters to be fitted");*/
      printf("LinFit: no parameters to be fitted\n"); exit(1);
   } 
   for (j=1;j<=mfit;j++) {             		/*Initialize the (symmetric) matrix.*/ 
      for (k=1;k<=mfit;k++) covar[j][k]=0.0; 
      beta[j][1]=0.0; 
   } 
   for (i=1;i<=ndat;i++) {             		/*Loop over data to accumulate coefficients of*/ 
      (*FouHarm)(x[i],afunc,x0,ma,Period); 	/*the normal equations.*/ 
      ym=y[i]; 
      if (mfit < ma) {                 		/*Subtract off dependences on known pieces*/ 
         for (j=1;j<=ma;j++)           		/*of the fitting function.*/ 
            if (!ia[j]) ym -= a[j]*afunc[j]; 
      } 
      sig2i=1.0/SQR(sig[i]); 
      for (j=0,l=1;l<=ma;l++) { 
         if (ia[l]) { 
            wt=afunc[l]*sig2i; 
            for (j++,k=0,m=1;m<=l;m++) 
               if (ia[m]) covar[j][++k] += wt*afunc[m]; 
            beta[j][1] += ym*wt; 
            b[j] = beta[j][1]; 
         } 
      } 
   } 
   for (j=2;j<=mfit;j++)               /*Fill in above the diagonal from symmetry.*/ 
      for (k=1;k<j;k++) 
         covar[k][j]=covar[j][k]; 
   ludcmp(covar,mfit,indx,d);          /*Matrix solution.*/
   lubksb(covar,mfit,indx,b); 
   for (j=0,l=1;l<=ma;l++) 
      if (ia[l]) a[l]=beta[++j][1];    /*Partition solution to appropriate coefficients a.*/ 
   *chisq=0.0; 
   for (i=1;i<=ndat;i++) {             /*Evaluate Chi2 of the fit.*/ 
      (*FouHarm)(x[i],afunc,x0,ma,Period);
      for (sum=0.0,j=1;j<=ma;j++) sum += b[j]*afunc[j];
      *chisq += SQR((y[i]-sum)/sig[i]); 
   } 
   covsrt(covar,ma,ia,mfit);           /*Sort covariance matrix to true order of fitting coefficients.*/
   free_vector(a,1,ma); free_vector(d,1,ma); free_ivector(indx,1,ma);  
   free_ivector(ia,1,ma); free_matrix(beta,1,ma,1,1); free_matrix(covar,1,ma,1,ma);  
} 
 
void covsrt(float **covar, int ma, int ia[], int mfit) 
/*Expand in storage the covariance matrix covar, so as to take into account parameters that are 
being held fixed. (For the latter, return zero covariances.)*/ 
/*Subroutine from Numerical Recipes in C*/ 
{ 
   float swap; 
   int i,j,k; 
 
   for (i=mfit+1;i<=ma;i++) 
      for (j=1;j<=i;j++) covar[i][j]=covar[j][i]=0.0; 
   k=mfit; 
   for (j=ma;j>=1;j--) { 
      if (ia[j]) { 
         for (i=1;i<=ma;i++) SWAP(covar[i][k],covar[i][j]) 
         for (i=1;i<=ma;i++) SWAP(covar[k][i],covar[j][i]) 
         k--; 
      } 
   } 
} 
 
void ludcmp(float **a, int n, int *indx, float *d) 
/*Given a matrix a[1..n][1..n], this routine replaces it by the LU decomposition of a rowwise 
permutation of itself. a and n are input. a is output, arranged as in equation (2.3.14) above; 
indx[1..n] is an output vector that records the row permutation effected by the partial 
pivoting; d is output as +/- 1 depending on whether the number of row interchanges was even 
or odd, respectively. This routine is used in combination with lubksb to solve linear equations 
or invert a matrix.*/ /*Subroutine from Numerical Recipes in C*/
{ 
   int i,imax,j,k; 
   float big,dum,sum,temp; 
   float *vv;                        	/*vv stores the implicit scaling of each row.*/ 
 
   vv=vector(1,n); 
   *d=1.0;                           	/*No row interchanges yet.*/ 
   for (i=1;i<=n;i++) {              	/*Loop over rows to get the implicit scaling*/ 
      big=0.0;                       	/*information*/ 
      for (j=1;j<=n;j++) 
         if ((temp=fabs(a[i][j])) > big) big=temp; 
      if (big == 0.0) {   	/*nrerror("Singular matrix in routine ludcmp");*/
         printf("Singular matrix in routine ludcmp\n"); exit(1);
      }                         	/*No nonzero largest element.*/
      vv[i]=1.0/big;                  	/*Save the scaling.*/ 
   } 
   for (j=1;j<=n;j++) {               	/*This is the loop over columns of Crout's method.*/ 
      for (i=1;i<j;i++) {             	/*This is equation (2.3.12) except for i = j.*/ 
         sum=a[i][j]; 
         for (k=1;k<i;k++) sum -= a[i][k]*a[k][j]; 
         a[i][j]=sum; 
      } 
      big=0.0;                         /*Initialize for the search for largest pivot element.*/ 
      for (i=j;i<=n;i++) {             /*This is i = j of equation (2.3.12) and i = j + 1...N*/ 
         sum=a[i][j];                  /*of equation (2.3.13).*/ 
         for (k=1;k<j;k++) 
            sum -= a[i][k]*a[k][j]; 
         a[i][j]=sum; 
         if ( (dum=vv[i]*fabs(sum)) >= big) { 
                        /*Is the figure of merit for the pivot better than the best so far?*/ 
            big=dum; imax=i; 
         } 
      } 
      if (j != imax) {               	/*Do we need to interchange rows?*/ 
         for (k=1;k<=n;k++) {        	/*Yes, do so...*/ 
            dum=a[imax][k]; a[imax][k]=a[j][k]; a[j][k]=dum; 
         } 
         *d = -(*d);                 	/*...and change the parity of d.*/ 
         vv[imax]=vv[j];             	/*Also interchange the scale factor.*/ 
      } 
      indx[j]=imax; 
      if (a[j][j] == 0.0) a[j][j]=TINY; 
      
      /*If the pivot element is zero the matrix is singular (at least to the precision of the 
      algorithm). For some applications on singular matrices, it is desirable to substitute 
      TINY for zero.*/ 
      
      if (j != n) {                  	/*Now, finally, divide by the pivot element.*/ 
         dum=1.0/(a[j][j]); 
         for (i=j+1;i<=n;i++) a[i][j] *= dum; 
      } 
   }                                 	/*Go back for the next column in the reduction.*/ 
   free_vector(vv,1,n); 
} 
 
void lubksb(float **a, int n, int *indx, float *b) 
/*Solves the set of n linear equations A . X = B. Here a[1..n][1..n] is input, not as the matrix 
A but rather as its LU decomposition, determined by the routine ludcmp. indx[1..n] is input 
as the permutation vector returned by ludcmp. b[1..n] is input as the right-hand side vector 
B, and returns with the solution vector X. a, n, and indx are not modified by this routine 
and can be left in place for successive calls with different right-hand sides b. This routine takes 
into account the possibility that b will begin with many zero elements, so it is efficient for use 
in matrix inversion.*/ /*Subroutine from Numerical Recipes in C*/
{ 
   int i,ii=0,ip,j; 
   float sum; 
 
   for (i=1;i<=n;i++) {             /*When ii is set to a positive value, it will become the index of*/ 
      ip=indx[i];					/*the first nonvanishing element of b. We now do the forward*/ 
      sum=b[ip];                    /*substitution, equation (2.3.6). The only new wrinkle is*/ 
      b[ip]=b[i];                   /* to unscramble the permutation as we go.*/ 
      if (ii) 
         for (j=ii;j<=i-1;j++) sum -= a[i][j]*b[j]; 
      else if (sum) ii=i;            	/*A nonzero element was encountered, so from now on we*/ 
      b[i]=sum;                      	/*will have to do the sums in the loop above.*/ 
   } 
   for (i=n;i>=1;i--) {              	/*Now we do the backsubstitution, equation (2.3.7).*/ 
      sum=b[i]; 
      for (j=i+1;j<=n;j++) sum -= a[i][j]*b[j]; 
      b[i]=sum/a[i][i];              	/*Store a component of the solution vector X.*/ 
   }                                 	/*All done!*/ 
} 

float SQR(float var)       /* Square function */
{ 
   return var*var; 
} 

/*All of the following subroutines are taken from Numerical Recipes in C*/
float *vector(long nl, long nh) 
/* allocate a float vector with subscript range v[nl..nh] */ 
{ 
   float *v; 
   v=(float *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(float))); 
   if (!v) {	/*nrerror("Allocation failure in vector()");*/
      printf("Allocation failure in vector()\n"); return 0;
   }
   return v-nl+NR_END; 
} 

int *ivector(long nl, long nh) 
/* allocate an int vector with subscript range v[nl..nh] */ 
{ 
   int *v; 
   v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int))); 
   if (!v) {	/*nrerror ("Allocation failure in ivector()");*/
      printf("Allocation failure in ivector()\n"); return 0;
   } 
   return v-nl+NR_END; 
} 

double *dvector(long nl, long nh) 
/* allocate a double vector with subscript range v[nl..nh] */ 
{ 
   double *v; 
   v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double))); 
   if (!v) {	/*nrerror("Allocation failure in dvector()");*/
      printf("Allocation failure in dvector()\n"); return 0;
   }
   return v-nl+NR_END; 
} 
 
float **matrix(long nrl, long nrh, long ncl, long nch)
/* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */ 
{ 
   long i, nrow=nrh-nrl+1,ncol=nch-ncl+1; 
   float **m; 
   m=(float **) malloc((size_t)((nrow+NR_END)*sizeof(float*)));  	/* allocate pointers to rows */ 
 
   if (!m) {	/*nrerror ("Allocation failure 1 in matrix()");*/
      printf("Allocation failure 1 in matrix()\n"); return 0;
   }
   m += NR_END; 
   m -= nrl; 
   m[nrl]=(float *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float)));	/* allocate rows and set pointers to them */ 
   if (!m[nrl]) {	/*nrerror ("Allocation failure 2 in matrix()");*/
         printf("Allocation failure 2 in matrix()\n"); return 0;
   }
   m[nrl] += NR_END; 
   m[nrl] -= ncl; 
   for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol; /* return pointer to array of pointers to rows */ 
   return m; 
} 

void free_cvector(unsigned char *v, long nl, long nh)
/* free an unsigned char vector allocated with cvector() */
{
   free((FREE_ARG) (v+nl-NR_END));
}
 
void free_vector(float *v, long nl, long nh) 
/* free a float vector allocated with vector() */ 
{ 
   free((FREE_ARG) (v+nl-NR_END)); 
} 

void free_ivector(int *v, long nl, long nh) 
/* free an int vector allocated with ivector() */ 
{ 
   free((FREE_ARG) (v+nl-NR_END)); 
} 

void free_dvector(double *v, long nl, long nh) 
/* free a double vector allocated with dvector() */ 
{ 
   free((FREE_ARG) (v+nl-NR_END)); 
} 
 
void free_matrix(float **m, long nrl, long nrh, long ncl, long nch) 
/* free a float matrix allocated by matrix() */ 
{ 
   free((FREE_ARG) (m[nrl]+ncl-NR_END)); 
   free((FREE_ARG) (m+nrl-NR_END)); 
}