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

#define DEG2RAD 0.01745329252
#define TWOPI 6.283185307
#define MAGIC_HEIGHT 220000.

/*  EVTH from cerfile. See CORSIKA manual for explanations */ 
typedef struct
{   char	EVTH[4];
    float	EvtNumber;
    float	PrimaryID;
    float	Etotal;   
    float	Thick0;   
    float	FirstTarget;
    float	zFirstInt;
    float	p[3];     
    float	Theta; 
    float	Phi; 
    float	NumRndSeq;
    float	RndData[10][3];
    float	RunNumber;
    float	DateRun;
    float	Corsika_version;
    float	NumObsLev;
    float	HeightLev[10]; 
    float	SlopeSpec;
    float	ELowLim;   
    float	EUppLim;   
    float	Ecutoffh;  
    float	Ecutoffm;  
    float	Ecutoffe;  
    float	Ecutoffg;  
    float	NFLAIN;
    float	NFLDIF;
    float	NFLPI0;
    float	NFLPIF;
    float	NFLCHE;
    float	NFRAGM; 
    float	Bx;
    float	By;
    float	EGS4yn;
    float	NKGyn;
    float	GHEISHAyn;
    float	VENUSyn;
    float	CERENKOVyn;
    float	NEUTRINOyn;
    float	HORIZONTyn;
    float	COMPUTER;
    float	ThetaMin;
    float	ThetaMax;
    float	PhiMin;
    float	PhiMax;
    float	CBunchSize;
    float	CDetInX,CDetInY;
    float	CSpacInX,CSpacInY;
    float	CLenInX,CLenInY; 
    float	COutput;
    float	AngleNorthX;
    float	MuonInfo;
    float	StepLength;
    float	CWaveLower;       
    float	CWaveUpper;       
    float	Multipl;       
    float	CorePos[2][20];   
    float       SIBYLL[2];
    float       QGSJET[2];
    float       DPMJET[2];
    float       VENUS_cross;
    float       mu_mult_scat;
    float       NKG_range;
    float       EFRCTHN[2];
    float       WMAX[2];
    float       rthin_rmax;
    float       viewcone_angles[2];
    float	dmmy2[119]; 
}   CerEventHeader;


/* Corsika Photon block: */
typedef struct
{   float w,       /*  photon wavelength (nm)      */
          x, y,    /*  (ground) impact point (cm)  */
          u, v,    /*  direction cosines           */
          t,       /*  arrival time (ns)           */
          h;       /*  production height (cm)      */
} photon;

/*****************************/

main(short argc, char *argv[])
{
float theta, phi;               // photons direction
double u, v, w;                  // direction cosines of photons
double source_x, source_y, source_z; // z is above sea level.
float xcenter = 0., ycenter = 0.; 
double cphz, norm;
float minradius, maxradius;     // area of photon impact points
long maxevents = 1;
long event, iphot;
FILE *out;

float dummy[272];
char runh[5], rune[5], evte[5];
CerEventHeader event_head;
photon cph;

float r, psi, maxpsi;

/*****************************/

if (argc < 4 || argc > 5)
  {
    printf ("\nUsage: cermaker source_x(cm) source_y(cm) source_z(cm) [events]\n\n");
    exit(-1);
  }

/* Source position with respect to telescope! */

sscanf (argv[1], "%lf", &source_x);
sscanf (argv[2], "%lf", &source_y);
sscanf (argv[3], "%lf", &source_z);

source_z += MAGIC_HEIGHT;

/*
sscanf (argv[4], "%f", &minradius);
sscanf (argv[5], "%f", &maxradius);
sscanf (argv[6], "%f", &maxpsi);
*/

minradius = 0.; maxradius = 900.; maxpsi = 360.;

if (argc == 5)
  sscanf (argv[4], "%ld", &maxevents);


memset ((char*)dummy,0,sizeof(dummy));
dummy[0] = 1;
dummy[2] = -1;
dummy[3] =  1;
dummy[4] =  220000;


memset ((char*)event_head.EVTH,0,sizeof(event_head));
event_head.EVTH[0] = 'E';
event_head.EVTH[1] = 'V';
event_head.EVTH[2] = 'T';
event_head.EVTH[3] = 'H';
event_head.PrimaryID = 1;
event_head.Etotal = 1.;
event_head.Theta = atan2(sqrt(source_x*source_x+source_y*source_y), source_z-MAGIC_HEIGHT);
event_head.Phi = atan2(-source_y,-source_x);
if (event_head.Phi < 0)
  event_head.Phi += TWOPI;
event_head.NumObsLev = 1; 
event_head.HeightLev[0] = 220000;
event_head.CERENKOVyn = 1;
event_head.Corsika_version = -1;
event_head.CBunchSize = 1;
event_head.CDetInX = 1;
event_head.CDetInY = 1;
event_head.COutput = 1;
event_head.CWaveLower = 290.;
event_head.CWaveUpper = 600.;
event_head.Multipl =1;
event_head.CorePos[0][0] = 0.;
event_head.CorePos[1][0] = 0.;

memset ((char*)&(cph.w),0,sizeof(cph));
cph.w = 500;
cph.t = 1000;

sprintf(evte, "EVTE");
sprintf(runh, "RUNH");
sprintf(rune, "RUNE");

out = fopen("cer000001","w");
if (!out)
  {
    printf ("\nERROR: Cannot open output file cer000001\n\n");
    exit(-1);
  }

fwrite (runh, strlen(runh), 1, out);
fwrite (dummy, sizeof(dummy), 1, out);

/* Here one event always means 39 photons! */

for (event = 0; event < maxevents; event++)
  {
    event_head.EvtNumber = event+1;
    fwrite (event_head.EVTH, sizeof(event_head), 1, out);

    for (iphot = 0; iphot < 39; iphot++)
      {
	psi = (double)rand()/(double)RAND_MAX*maxpsi*DEG2RAD;
	r  = sqrt( (double)rand()/(double)RAND_MAX*
	  (maxradius*maxradius-minradius*minradius) + 
	    minradius*minradius);

	cph.x = xcenter+r*cos(psi);
	cph.y = ycenter+r*sin(psi);

	cphz = MAGIC_HEIGHT;        /* Above sea level */


	/* Photon direction: */
	u = (double)cph.x - source_x;
	v = (double)cph.y - source_y;
	w =         cphz  - source_z;

	norm = sqrt(u*u+v*v+w*w);

	cph.u = (float)(u/norm);
	cph.v = (float)(v/norm);
        w = (float)(w/norm);

	cph.h =  source_z;

	fwrite(&(cph.w),sizeof(cph),1,out);
      }

    fwrite (evte, strlen(evte), 1, out);
    fwrite (dummy, sizeof(dummy), 1, out);
  }

fwrite (rune, strlen(rune), 1, out);
fwrite (dummy, sizeof(dummy), 1, out);
    
fclose(out);

return;

}