/* ----------------------------------------------------------------------
 *
 *  Created:  Thu May  7 16:24:22 1998
 *  Author:   Jose Carlos Gonzalez
 *  Purpose:  Program for reflector simulation
 *  Notes:    See files README for details
 *
 * ----------------------------------------------------------------------
 *
 *  Modified: Tue Jan 16 10:10:10 2001
 *  Authors:  Denis Bastieri + Ciro Bigongiari
 *  Purpose:  Program for reflector simulation (Single file input)
 *
 * ----------------------------------------------------------------------
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "version.h"
#include "diag.h"
#include "init.h"
#include "header.h"

FILE *rflfile = NULL,		/*  reflector (output) file	*/
     *cerfile = NULL;		/*  current cerfile		*/

cphoton *CPhotons = NULL;		  /*  cphoton array	*/
static photon Photons[PH_IN_DATABLOCK];   /*  photon array	*/

static long readp = 0L;		/*  read ptr (on cerfile)	*/
static long writep = 0L;	/*  write ptr (on rflfile)	*/


extern CerHeader *cheadp;	/*  var inited in header.c	*/
extern RflHeader *rheadp;	/*  var inited in header.c	*/

/*  Prototypes  */
void clean(void);
FILE *GetNextFile(char *cername);
static int GetEvent(void);
static int ProcessEvent(CerHeader *cheadp, FILE *cerfile, FILE *rflfile);

void main(void)
{   long event = 0L;		/*  event counter	*/

    /*  Read init & geometry parms, init files and vars */
    init(NULL);

    /*	Processing loop	*/
    while(event < max_Events && GetEvent())
    {	if (ProcessEvent(cheadp, cerfile, rflfile)) event++;
	if (event % SHOW_ME == 0) Log(INFO_EVENT_LOG,
	    cheadp->RunNumber, cheadp->EvtNumber, cheadp->Etotal);  }

    /*  Writing final flags  */
    fwrite(FLAG_END_OF_FILE, SIZE_OF_FLAGS, 1, rflfile);

    /*  Close file  */
    Log(RFLF_CLOSE_LOG);
    fclose(rflfile);

    /*  Clean memory and exit  */
    clean();
    Message(RFL__EXIT__MSG, QUOTE(PROGRAM), QUOTE(VERSION));

}   /*  end of main  */

static int ProcessEvent(CerHeader *cheadp, FILE *cerfile, FILE *rflfile)
{   extern int absorption(float wlen, float height, float theta);
    extern int ph2cph(photon *ph, cphoton *cph);
    extern float Omega[3][3];
    extern float OmegaI[3][3];
    extern float OmegaCT[3][3];
    extern float OmegaICT[3][3];
    extern void makeOmega(float theta, float phi);
    extern void makeOmegaI(float theta, float phi);

    /*  Various counters: phs = absphs + refphs[0..3] + cphs  */
    long phs,		/*  Number of incoming photons	*/
	 absphs,	/*  Photons absorbed		*/
	 refphs[4],	/*  Photons not reflected	*/
	 cphs;		/*  Number of cphotons		*/

    FILE *tmpf=NULL;	/*  Temp fp to handle o/f	*/
    size_t read;	/*  items read: != 1 on error	*/
    int overflow=0,	/*  how many o/f on cphs	*/
	ref_type,	/*  ret value from reflection	*/
	ph;		/*  photon number		*/

    float first = 1e8f,	/*  Photon arrival times	*/
	  last  = 0;

    /*	photon quantities  */
    float wlen,		/*  photon wavelength	*/
	  theta;	/*  photon zenith angle	*/

    /*	Reset counters and set fileptr  */
    phs= absphs= refphs[0]= refphs[1]= refphs[2]= refphs[3]= cphs= 0;
    writep = ftell(rflfile);

    /*  Translate the EVTH from cerfile to rflfile  */
    TranslateHeader(rheadp, cheadp);

    /*  Calculate OmegaCT matrices  */
    if (is_Fixed_Target)
    {	makeOmega(fixed_Theta, fixed_Phi);
	makeOmegaI(fixed_Theta, fixed_Phi);   }
    else
    {	makeOmega(cheadp->Theta, cheadp->Phi);
	makeOmegaI(cheadp->Theta, cheadp->Phi);   }
    memcpy( OmegaCT, Omega, 9*sizeof(float) );
    memcpy( OmegaICT, OmegaI, 9*sizeof(float) );

    /*  Loop on data blocks  */
    while(1 == (read = fread(Photons, sizeof(Photons), 1, cerfile))
      &&  strncmp((char *)Photons, "EVTE", 4))

    /*  Loop on phs inside block: exit when wlength is 0  */
    {	for (ph=0; ph<PH_IN_DATABLOCK; ph++)
	{   if (Photons[ph].w == 0) break;

	    CPhotons[cphs].w = Photons[ph].w;
	    Photons[ph].w = wlen = (float) fmod(Photons[ph].w, 1000.);
	    Photons[ph].x -= cheadp->CorePos[0][0];
	    Photons[ph].y -= cheadp->CorePos[1][0];

	    /*  Increment number of photons	*/
	    phs++;
	    
	    /*  Calculate some quantities  */
	    theta = (float) acos(sqrt(
		    1.f - Photons[ph].u*Photons[ph].u
			- Photons[ph].v*Photons[ph].v));
	    /*	Check absorption  */

	    if (absorption(wlen, Photons[ph].h, theta))
		absphs++;

	    /*	Check reflection  */
	    else if (0 != (ref_type =
			   ph2cph(&Photons[ph], &CPhotons[cphs]))) 
		refphs[ref_type-1]++; 

	    else    /*	Photon passed	*/
	    {
		/*  Update first/last arrival times  */
		if (first > CPhotons[cphs].t) first = CPhotons[cphs].t;
		if ( last < CPhotons[cphs].t)  last = CPhotons[cphs].t;

		/*  Update cphs	*/
		if (++cphs == NR_OF_CPHOTONS)	/*  Overflow  */
		{
		    /*  if it is the first o/f open a tempfile  */
		    if (overflow++ == 0)
		    {	Log("Spooling... ");
			tmpf = tmpfile();
			if (tmpf == NULL) FatalError(TEMP_ERROR_FTL);  }

		    /*  Write now the cphoton array	*/
		    fwrite(CPhotons, sizeof(cphoton), NR_OF_CPHOTONS, tmpf);

		    /*  Reset the cphs counter	*/
		    cphs = 0;

		}   /*  if overflow  */
	    }  /*  else (=Photon passed)  */
	}   /*  end of loop inside datablock	*/
    }   /*  end of loop on datablocks		*/

    /*  Check if there was an error while reading cerfile  */
    if (read != 1) fseek(rflfile, writep, SEEK_SET);

    else    /*  no error: write the new event	*/

    {	/*  Write "start of event" flag  */
	fwrite(FLAG_START_OF_EVENT, SIZE_OF_FLAGS, 1, rflfile);

	/*  Update arrival times  */
	rheadp->TimeFirst = first;
	rheadp->TimeLast  = last;

	/*  Update RflHeader with info on ph/cph nrs and write it  */
	rheadp->CORSIKAPhs   = phs;
	rheadp->AtmAbsPhs    = absphs;
	rheadp->MirrAbsPhs   = refphs[0];
	rheadp->OutOfMirrPhs = refphs[1];
	rheadp->BlackSpotPhs = refphs[2];
	rheadp->OutOfChamPhs = refphs[3];
	rheadp->CPhotons     = (long) overflow * NR_OF_CPHOTONS + cphs;
	fwrite(rheadp, sizeof(RflHeader), 1, rflfile);

	/*  If there was an overflow, append data from tempfile   */
	if (overflow)
	{
	    /*  Unload data from CPhotons  */
	    fwrite(CPhotons, sizeof(cphoton), cphs, tmpf);

	    /*  Transfer data  */
	    fseek(tmpf, 0L, SEEK_SET);	/*  Start from the beginning  */
	    while (overflow--)
	    {	fread (CPhotons, sizeof(cphoton), NR_OF_CPHOTONS, tmpf);
		fwrite(CPhotons, sizeof(cphoton), NR_OF_CPHOTONS, rflfile); }

	    /*  Reload data in CPhotons    */
	    fread (CPhotons, sizeof(cphoton), cphs, tmpf);

	    /*	Close (and delete) temp file  */
	    fclose(tmpf);  }

	/*  Write (remaining) cphotons  */
	fwrite(CPhotons, sizeof(cphoton), cphs, rflfile);

	/*  Write "end of event" flag  */
	fwrite(FLAG_END_OF_EVENT, SIZE_OF_FLAGS, 1, rflfile);

	Debug(">>> %6ld = %6ld + %6ld + %6ld + %6ld + %6ld + %6f\n",
	      phs, absphs,
	      refphs[0], refphs[1], refphs[2], refphs[3],
	      rheadp->CPhotons); }

    return read == 1;
}   /*	end of ProcessEvent  */

static int GetEvent(void)
{   int found = FALSE,		/*  event found		*/
	isWrong = FALSE;	/*  cerfile is wrong	*/
    static int newFile = TRUE;  /*  if TRUE, check if cerfile is valid	*/

    do
    {	/*  In case the just-opened file is a valid cerfile,
	    starting with a RUNH, loop until a valid event is found:
	    id est with: first_Event <= EvtNumber <= last_Event
	    and             low_Ecut <=  Etotal   <= high_Ecut
	    If the search was successful, "found" is set to TRUE.
	    If there are reading errors, "isWrong" is set to TRUE. */

	if (newFile
	 && (1 != fread(cheadp, sizeof(CerHeader), 1, cerfile)
	     || strncmp(cheadp->EVTH, "RUNH", 4)))
	{   isWrong = TRUE;   }

	else do
	{   /*  Push fileptr in case it is a "EVTH"  */
	    readp = ftell(cerfile);

	    /*  Error: exit loop  */
	    if (1 != fread(cheadp, sizeof(CerHeader), 1, cerfile))
	    {	isWrong = TRUE;    break;   }

	    /*  Ok: set found at TRUE and exit loop  */
	    if (strncmp(cheadp->EVTH, "EVTH", 4) == 0
	     && first_Event <= cheadp->EvtNumber
	     &&		       cheadp->EvtNumber <= last_Event
	     &&    low_Ecut <= cheadp->Etotal
	     &&		       cheadp->Etotal <= high_Ecut)
	    {	found = TRUE;   }

	    /*  File is finished: exit loop  */
	    else if (strncmp(cheadp->EVTH, "RUNE", 4) == 0)
		break;

	} while(found == FALSE);

	/*  If there was an error, log it  */
	if (isWrong) Log(CERF_ERROR_LOG, cername);

	/*  If a new event was not found, get, if any, a new cerfile.
	    "newFile" is set in this case, to check the validity of
	    just-opened cerfile.  If GetNextFile return a NULL (no
	    new cerfile), then exit thoroughly.  */

	if (found) newFile = FALSE;
	else
	{   if ((cerfile=GetNextFile(cername)) == NULL) break;
	    newFile = TRUE;   }
	
    } while(found == FALSE);

    return found;
}   /*  end of GetEvent  */

/*  Files are all CER-like files (cerfiles).
    Filenames are written one per line and may be followed
    by a ",FIRST_EVT,LAST_EVT" directive, such as:
 -> cer19999 15 167
    telling to deal with events from 15 to 167 in file "cer19999"
    The file list may be appended to the input file or can be stored
    in a separate file pointed by an '@' directive contained in the
    input file.  In both cases, file referral follows the tag
    "cer_files".  Files not accessible are thoroughly skipped.
    GetNextFile gets and opens the next readable cerfile.  */

FILE *GetNextFile(char *cername)
{   FILE *f = NULL;		/*  return value (cerfile ptr)	*/
    char *value_ptr;		/*  ptr at parm value	*/
    extern char line[];		/*  white chars (init)	*/
    extern char whites[];	/*  parsing buf. (init)	*/

    /*  Close, if any, the previous cerfile, writing "END_OF_RUN".  */
    if (cerfile)
    {	fclose(cerfile);
	fwrite(FLAG_END_OF_RUN, SIZE_OF_FLAGS, 1, rflfile);  }

    /*  Parse next line in filelist  */
    do
    {	if (fgets(line, LINE_MAX_LENGTH, filelist) == NULL) break;
	if ((value_ptr=strtok(line, whites)) == NULL) continue;

	/*  If you found a line with some meaning, try to open the file  */
	if ((f=fopen(value_ptr, "r")) == NULL)
	    Message(CERF_NFND__MSG, value_ptr);

	else	    /*  Cerfile is readable  */
	{   /*  Store its name in "cername"  */
	    Log(CERF_OPEN__LOG, strcpy(cername, value_ptr));

	    /*  Read, if any, event bounds   */
	    if ((value_ptr=strtok(NULL, whites)) == NULL)
	    {	first_Event = 0;
		last_Event  = 1000000;   }
	    else
	    {	first_Event = atol(value_ptr);
		value_ptr = strtok(NULL, whites);
		last_Event = value_ptr ? atol(value_ptr) : 1000000;  }

	    /*  Check boundaries  */
	    if (first_Event > last_Event)
	    {	Error(EVTN_WRONG_ERR, first_Event, last_Event, cername);
		first_Event = 0;
		last_Event  = 1000000;  }}}

    while (f == NULL);	/*  Loop until a readable file is found  */

    /*	If a new cerfile is met, write "START_OF_RUN"  */
    if (f) fwrite(FLAG_START_OF_RUN, SIZE_OF_FLAGS, 1, rflfile);

    return f;
}   /*  end of GetNextFile  */

/*  This function frees up allocated memory before exiting  */

void clean(void)
{   Message(MEM__FREE__MSG);

    if (ct_data)
    {	free(ct_data);		Log(VECT_FREE__LOG, "ct_data");  }

    if (Reflectivity[0])
    {	free(Reflectivity[0]);	Log(VECT_FREE__LOG, "Reflectivity[0]");  }
    if (Reflectivity[1])
    {	free(Reflectivity[1]);	Log(VECT_FREE__LOG, "Reflectivity[1]");  }

    if (AxisDeviation[0])
    {	free(AxisDeviation[0]);	Log(VECT_FREE__LOG, "AxisDeviation[0]");  }
    if (AxisDeviation[1])
    {	free(AxisDeviation[1]);	Log(VECT_FREE__LOG, "AxisDeviation[1]");  }

    if (ct_Focal)
    {	free(ct_Focal);		Log(VECT_FREE__LOG, "ct_Focal");  }

    if (CPhotons)
    {	free(CPhotons);		Log(VECT_FREE__LOG, "CPhotons");  }
}   /*  end of clean  */
