/* compiles by `gcc ucr.c -I. -Lrdmc -lrdmc -lm -o ucr -O2 -Wall -static' */

#include <rdmc/rdmc.h>
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>

#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>

#define VERSION "0.05"
#define MAX_IN 256
#define MAX_STEPS 20
#define PRECISION 1e-8
#define PI 3.141592653589793
#define SPEEDOFLIGHT 0.299792458

char randomize=0;
char curved=0;
char rphi=0;
char rmpri=0;
char rmusr=0;
char oms=0;
char ohm=0;
char msn=0;
char remr=0;
char remrflag=0;
char rdmcF2k=0;

int RUNNUM=0;
int TRNUM=1;
int OVERSAMPLE=1;
int insert=0;

double EARTHR=6.4e6;
double LENGTH=800;
double RADIUS=400;
double DEPTH=1730;
double HEIGHT=2834;
double DCORR=0;
double REMR=0;
double cutth=0;
double cutfe=0;
double llow=0;
double trigw=-1;
double avtime=0;
double alow=0.212/1.2;
double blow=0.251e-3/1.2;

double etrans=0.;
double dtrans=0.;
double gmaxdtrans=2000.;
double transheight=25000.;

struct stat outstat;
char outfile[MAX_IN]="../output/outfile.gz";


double alpha(double z)
{
  int i;
  double x, dx, x1, x2, f1, f2, df, f;

  x1=0; x2=PI/2;
  f1=x1+0.5*sin(2*x1)-z;
  f2=x2+0.5*sin(2*x2)-z;
  x=x1;

  if(f1*f2>0){
    printf("Warning: Root is not bracketed\n");
    return x;
  }

  for(i=0; i<MAX_STEPS; i++){
    f=x+0.5*sin(2*x)-z;
    if(f==0) return x;
    else if(f1*f<0){ x2=x; f2=f; }
    else{ x1=x; f1=f; }
    df=1+cos(2*x);
    if(df!=0) dx=-f/df; else dx=0;
    if(x+dx<x1 || x+dx>x2 || df==0) dx=(x2-x1)/2;
    x+=dx;
    if(fabs(dx)<PRECISION) break;
  }
  if(i==MAX_STEPS) printf("Warning: Failed to achieve accuracy < %g in %i steps\n", PRECISION, i);

  return x;
}


void delta(double *dx, double *dy, double costh)
{
  double x, y, a, b, l, r, sinth, w;

  *dx=0, *dy=0;
  a=rand(); b=rand();
  a/=RAND_MAX, b/=RAND_MAX;
  l=LENGTH, r=RADIUS;

  sinth=sqrt(1-costh*costh);
  w=costh+sinth*(2*l)/(PI*r);

  if(w*a < 0.5*costh){
    x=-r*sin(alpha(PI*(0.5-a*w/costh)));
    y=(2*b-1)*sqrt(r*r-x*x);
    x=x*costh-(l/2)*sinth;
  }
  else if(w*a <= w-0.5*costh){
    x=r*(PI/2)*(a*w-0.5*costh)-(l/2)*sinth;
    y=(2*b-1)*r;
  }
  else {
    x=r*sin(alpha(PI*(0.5-(1-a)*w/costh)));
    y=(2*b-1)*sqrt(r*r-x*x);
    x=x*costh+(l/2)*sinth;
  }

  *dx=x; *dy=y;

  return;
}


int process(char *filename, mcfile *fo)
{
  long int i, j;
  mcfile *fi;
  static array ar, ar_new;  // static needed since not using rdmc_init_array
  static mevt ev, ev_save, ev_new;   // or rdmc_init_mevt
  int flag, oversample, error=0, saved_i;
  double deltax, deltay, deltaz, dx, dy, dz, dphi=0, aphi, e;
  double phi, cosp, sinp, costh, sinth, sina, cosa, xc, yc, zc, b, c, r;
  double axx, axy, axz, ayx, ayy, ayz, azx, azy, azz, rx, ry, rz, px, py, pz;
  double maxdtrans=0, mindtrans=0, mintime=0;
  double ptdist=0, ptdistx=0, ptdisty=0, holdtheta=0;
  static int num=0;
  static long int enr=0, evtri=0, numEM=0, numTR=0, newEM=0, newTR=0;
  static double evtime=0, uttime=0;

  if((fi=rdmc_mcopen(filename, "r", F2000_F))==NULL){
    printf("cannot open input file `%s', exiting\n", filename);
    return 1;
  }

  if(!num) rdmc_mccpfp(fi, fo);
  rdmc_rarr(fi, &ar);
  if(!num){
    if(stat(outfile, &outstat)) exit(1);
    rdmc_cp_array(&ar_new, &ar);
    if(rmusr) while(ar_new.n_user>0) rdmc_del_user_def(&ar_new, 0);

    if(randomize && ar_new.n_user==2){
      static array_hdef_t nu;     // or rdmc_init_array_hdef
      nu.id=2;
      strcpy(nu.tag, "RANDXY");
      nu.nwords=3;
      strcpy(nu.words[0], "COREX");
      strcpy(nu.words[1], "COREY");
      strcpy(nu.words[2], "RPHI");
      rdmc_add_user_def(&ar_new, &nu, 2);
    }
    rdmc_warr(fo, &ar_new);
    num=1;
  }

  while(1){
    error=rdmc_revt(fi, &ar, &ev);
    if(error==EOF){
      printf("EOF numEM=%li numTR=%li newEM=%li newTR=%li\n", numEM, numTR, newEM, newTR);
      break;
    }
    if(error==RDMC_ILF){
      printf("File format error\n");
      break;
    }
    if(error){
      printf("Format error %i\n", error);
      break;
    }

    if(RUNNUM) ev.nrun=RUNNUM;

    numEM++; numTR+=ev.ntrack;
    if(insert){
      // insert transverse momentum muon, but only if there's already muons, also get min time of all muons
      mintime=1e9;
      for(i=0;i<ev.ntrack;i++){
	if(ev.gen[i].id==5||ev.gen[i].id==6){  // grabs muons
	  if(ev.gen[i].t<mintime)mintime=ev.gen[i].t;
	}
      }
      if(mintime==1e9) continue;

      // Get random separation
      ptdistx=ptdisty=ptdist=0.;
      maxdtrans=gmaxdtrans;
      mindtrans=-dtrans;
      if(dtrans>0){          // Positive dtrans is a fixed separation
	maxdtrans=1.1*dtrans;
	mindtrans=0.9*dtrans;
      }
      while((ptdist<mindtrans)||(ptdist>maxdtrans)){
	ptdist=((float)rand()/(float)RAND_MAX)*(maxdtrans-mindtrans)+mindtrans;
	holdtheta=2.*3.141592654*((float)rand()/(float)RAND_MAX);
	ptdistx=ptdist*cos(holdtheta);
	ptdisty=ptdist*sin(holdtheta);
      }

      static mtrack hold;    // or rdmc_init_mtrack
      hold.x=ptdistx;
      hold.y=ptdisty;
      hold.z=ev.gen[0].z;
      hold.t=mintime;
      hold.id=5;
      hold.tag=ev.ntrack+1;
      hold.costh=ev.gen[0].costh;
      hold.phi=ev.gen[0].phi;
      hold.e=etrans*1.e3;    // track energy is in MeV, but it's input is in GeV
      rdmc_add_gen(&ev, &hold, ev.ntrack);
      numTR++;
    }

    if(trigw>=0){
      r=rand(); r=-log(r/RAND_MAX)*avtime;
      evtime+=r; uttime+=r;
    }

    if(oms || msn){
      flag=0;
      if(ohm){
	e=0;
	saved_i=-1;
      }
      for(i=1; i<ev.ntrack; i++){
	if((ev.gen[i].id==5) || (ev.gen[i].id==6) ||
	   (msn && (ev.gen[i].id==4 || ev.gen[i].id==201 || ev.gen[i].id==202 || ev.gen[i].id==203 ||
		  ev.gen[i].id==204 || ev.gen[i].id==205 || ev.gen[i].id==206))){
	  if(ohm){
	    if(ev.gen[i].e>e){
	      if(saved_i!=-1){
		rdmc_del_gen(&ev, saved_i);
		i--;
	      }
	      e=ev.gen[i].e;
	      saved_i=i;
	    }
	    else {
	      rdmc_del_gen(&ev, i);
	      i--;
	    }
	  }
	  flag=1;
	}
	else {
	  rdmc_del_gen(&ev, i);
	  i--;
	}
      }
      if(!flag) continue;
    }

    if(llow!=0){
      for(i=1; i<ev.ntrack; i++){
	if((ev.gen[i].id==5) || (ev.gen[i].id==6)){
	  flag=0;
	  cutfe=ev.gen[i].costh;
	  if(cutfe==0) flag=1;
	  else{
	    cutfe=blow*llow/cutfe;
	    if(cutfe>100) cutfe=100;
	    cutfe=(exp(cutfe)-1)*alow/blow;
	    if(cutfe>ev.gen[i].e/1.e3) flag=1;
	  }
	  if(flag){
	    rdmc_del_gen(&ev, i);
	    i--;
	  }
	}
      }
    }

    if(rmpri) if(ev.ntrack>0) rdmc_del_gen(&ev, 0);
    if(rmusr) while(ev.nuser>0) rdmc_del_user(&ev, 0);

    if(ev.ntrack<TRNUM){
      continue;
    }

    if(rdmcF2k){
      for(i=0; i<ev.ntrack; i++){ if(ev.gen[i].id==14) ev.gen[i].id=3001; ev.gen[i].tag=i; }
      for(i=1; i<ev.ntrack; i++) ev.gen[i].parent=0;
    }

    if(remr){
      deltax=ev.gen[0].x;
      deltay=ev.gen[0].y;
      if(!remrflag){
	if(ev.ntrack>1) deltaz=ev.gen[1].z;
	else deltaz=DEPTH;
      }
      else deltaz=REMR;
      for(j=0; j<ev.ntrack; j++){
	ev.gen[j].z-=deltaz-HEIGHT;
	ev.gen[j].x-=deltax;
	ev.gen[j].y-=deltay;
      }
    }

    oversample=OVERSAMPLE;
    if(curved==4) if(ev.gen[0].costh<=cutth) oversample*=2;
    if(oversample!=1) rdmc_cp_mevt(&ev_save, &ev);
    for(i=0; i<oversample; i++){
      if(i>0) rdmc_cp_mevt(&ev, &ev_save);
      enr++;
      ev.enr=enr;

      if(rphi){
	dphi=rand(); dphi/=RAND_MAX; dphi*=2*PI;
	for(j=0; j<ev.ntrack; j++){
	  aphi=ev.gen[j].phi+dphi;
	  if(aphi>2*PI) aphi-=2*PI;
	  ev.gen[j].phi=aphi;
	  dx=ev.gen[j].x*cos(dphi)-ev.gen[j].y*sin(dphi);
	  dy=ev.gen[j].x*sin(dphi)+ev.gen[j].y*cos(dphi);
	  ev.gen[j].x=dx;
	  ev.gen[j].y=dy;
	}
      }

      if(randomize){
	for(j=0; j<ev.ntrack; j++){
	  ev.gen[j].z-=HEIGHT-DEPTH;    // Depth correction: surface -> AMANDA coordinates
	}

	phi=ev.gen[0].phi; costh=ev.gen[0].costh;

	if(!curved){
	  delta(&dx, &dy, costh); dx/=costh;
	  dx+=DEPTH*sqrt(1/(costh*costh)-1);
	  deltax=cos(phi)*dx-sin(phi)*dy;
	  deltay=sin(phi)*dx+cos(phi)*dy;
	}
	else{
	  sinth=sqrt(1-costh*costh);
	  sinp=sinth/(1-DEPTH/EARTHR);  // angle in the detector frame
	  if(sinp>1+PRECISION){
	    rdmc_mcclose(fi);
	    printf("ERROR: input file `%s' is incompatible with the curved option\n\
Please set the SCURV option to true in the INPUTS file\n", filename);
	    return 1;
	  }
	  else if(sinp>1) sinp=1;
	  cosp=sqrt(1-sinp*sinp);

	  delta(&dx, &dy, cosp); dx/=costh;
	  deltax=cos(phi)*dx-sin(phi)*dy;
	  deltay=sin(phi)*dx+cos(phi)*dy;

	  switch(curved){
	  case 1: break;
	  case 2: cosp=-cosp; break;
	  case 3: if(rand()<RAND_MAX/2) cosp=-cosp; break;
	  case 4: if(costh<=cutth) if(i%2==1) cosp=-cosp; break;
	  default: break;
	  }

	  sina=sinp*costh-cosp*sinth; if(sina>1) sina=1; // angle difference between the two frames
	  cosa=cosp*costh+sinp*sinth; if(cosa>1) cosa=1; else if(cosa<-1) cosa=-1;

	  cosp=cos(phi);
	  sinp=sin(phi);

	  axx=cosp*cosp*(cosa-1)+1;     // transformation matrix computation
	  axy=cosp*sinp*(cosa-1);
	  axz=cosp*sina;
	  ayx=axy;
	  ayy=sinp*sinp*(cosa-1)+1;
	  ayz=sinp*sina;
	  azx=-axz;
	  azy=-ayz;
	  azz=cosa;

	  rx=EARTHR*sina*cosp;          // shift vector computation
	  ry=EARTHR*sina*sinp;
	  rz=EARTHR*(cosa-1);
	}

	if(randomize && ar.n_user==3){
	  static mevt_special_t nu;     // or rdmc_init_mevt_special
	  nu.id=2;
	  nu.nval=3;
	  nu.val[0]=deltax;
	  nu.val[1]=deltay;
	  nu.val[2]=dphi;
	  rdmc_add_user(&ev, &nu, 2);
	}

	for(j=0; j<ev.ntrack; j++){
	  if(!curved){
	    ev.gen[j].x+=deltax;
	    ev.gen[j].y+=deltay;
	    ev.gen[j].z-=DCORR;
	  }
	  else{
	    dx=ev.gen[j].x+deltax;      // coordinates in the local CORSIKA frame
	    dy=ev.gen[j].y+deltay;
	    dz=ev.gen[j].z;

	    xc=axx*dx+axy*dy+rx;        // coordinates on the tangent plane in the detector frame
	    yc=ayx*dx+ayy*dy+ry;
	    zc=azx*dx+azy*dy+rz+dz;     // z is additionally corrected by HEIGHT-DEPTH

	    phi=ev.gen[j].phi; costh=ev.gen[j].costh;
	    sinth=sqrt(1-costh*costh);
	    cosp=cos(phi); sinp=sin(phi);

	    dx=sinth*cosp;              // direction in the CORSIKA frame
	    dy=sinth*sinp;
	    dz=costh;

	    px=axx*dx+axy*dy+axz*dz;    // direction in the detector frame
	    py=ayx*dx+ayy*dy+ayz*dz;
	    pz=azx*dx+azy*dy+azz*dz;

	    costh=pz; if(costh>1) costh=1; else if(costh<-1) costh=-1;
	    sinth=sqrt(1-costh*costh);
	    if(sinth!=0){
	      sinp=py/sinth;
	      cosp=px/sinth; if(cosp>1) cosp=1; else if(cosp<-1) cosp=-1;
	    }
	    phi=acos(cosp); if(sinp<0) phi=2*PI-phi; if(phi>=2*PI) phi-=2*PI;

	    ev.gen[j].costh=costh;
	    ev.gen[j].phi=phi;

	    b=xc*px+yc*py+(zc+EARTHR-DEPTH)*pz;
	    c=xc*xc+yc*yc+(zc+EARTHR-DEPTH)*(zc+EARTHR-DEPTH)-(EARTHR-DCORR)*(EARTHR-DCORR);
	    r=b-sqrt(b*b-c);            // shift along (px,py,pz) to reach the surface from the tangent plane

	    ev.gen[j].x=xc-px*r;        // final coordinates in the detector frame
	    ev.gen[j].y=yc-py*r;
	    ev.gen[j].z=zc-pz*r;
	    ev.gen[j].t+=r/SPEEDOFLIGHT;
	    if(ev.gen[j].length!=-1) ev.gen[j].length+=r;

	  }
	}
      }

      if(trigw>=0) if(ev.ntrack>0){
	phi=ev.gen[0].phi; costh=ev.gen[0].costh;
	dx=ev.gen[0].x; dy=ev.gen[0].y; dz=ev.gen[0].z;
	sinp=sin(phi);
	cosp=cos(phi);
	sinth=sqrt(1-costh*costh);
	r=dx*sinth*cosp+dy*sinth*sinp+dz*costh;
	for(j=0; j<ev.ntrack; j++) ev.gen[j].t-=r/SPEEDOFLIGHT;
	if(uttime>=trigw*1.e-6){
	  if(ev_new.ntrack>0) newEM++; newTR+=ev_new.ntrack;
	  if(stat(outfile, &outstat)) exit(1);
	  rdmc_wevt(fo, &ev_new, &ar_new);
	  evtri=0; uttime=0;
	}
	if(evtri==0){
	  j=floor(evtime/(24*60*60));
	  ev.secs=floor(evtime-j*(24*60*60));
	  ev.nsecs=floor(1.e9*(evtime-j*(24*60*60)-ev.secs));
	  ev.mjd+=j;
	  rdmc_cp_mevt(&ev_new, &ev);
	}
	else{
	  for(j=0; j<ev.ntrack; j++){
	    ev.gen[j].t+=uttime*1.e9;
	    ev.gen[j].tag+=evtri;
	    if(ev.gen[j].parent>=0) ev.gen[j].parent+=evtri;
	    rdmc_add_gen(&ev_new, &ev.gen[j], evtri+j);
	  }
	}
	evtri+=ev.ntrack;
	continue;
      }

      newEM++; newTR+=ev.ntrack;
      if(stat(outfile, &outstat)) exit(1);
      rdmc_wevt(fo, &ev, &ar_new);
    }
  }

  rdmc_mcclose(fi);
  if(error==EOF) error=0;
  else if(error==RDMC_ILF) error=1;
  else if(error) error=1;
  return error;
}


int my_input(char *s, FILE *fc)
{
  int i;

  s[0]=' ';
  for(i=1; i<MAX_IN-1; i++){
    if(fscanf(fc, "%c", &s[i])!=1){
      s[i]='\n', s[i+1]='\0';
      return 1;
    }
    if(s[i]=='\n'){
      i++;
      break;
    }
  }
  s[i]='\0';
  return 0;
}

int main(int arg_c, char *arg_a[])
{
  FILE *fc;
  mcfile *fi, *fo;
  static array ar;    // all static variables initialized with 0
  static mevt ev;     //   this ensures all pointers are NULL
  char setcmt=0;
  char s[MAX_IN];
  char cmtfile[MAX_IN];
  int i=0, flag=0, length, error=0, errf=0, j, jmax;
  long int test_num=0, evts=0;
  double AREASUM=0, FLUXSUM=0, SHOWERS=0, LIFETIME=0;
  double test, test_cos=0, test_phi, test_x, test_y, dx, dy;

  if(arg_c==1){
    printf("This program processes CORSIKA f2k files. Use `%s -h' for help\n", arg_a[0]);
    return 1;
  }

  for(i=1; i<arg_c; i++){
    if(*arg_a[i]=='-'){

      if(!strncmp(arg_a[i], "-out=", 5)){
	length=strlen(arg_a[i]+5);
	if(length>=MAX_IN){
	  printf("File name too long\n");
	  return 1;
	}
	if(length<=0){
	  printf("File name too short\n");
	  return 1;
	}
	strcpy(outfile, arg_a[i]+5);
	continue;
      }

      if(!strncmp(arg_a[i], "-tr=", 4)){
	TRNUM=atoi(arg_a[i]+4);
	if(TRNUM<1) TRNUM=1;
	continue;
      }

      if(!strncmp(arg_a[i], "-over=", 6)){
	OVERSAMPLE=atoi(arg_a[i]+6);
	if(OVERSAMPLE<1) OVERSAMPLE=1;
	randomize=1;
	continue;
      }

      if(!(strcmp(arg_a[i], "-r") && strcmp(arg_a[i], "-randomize"))){
	randomize=1;
	continue;
      }

      if(!(strcmp(arg_a[i], "-c") && strcmp(arg_a[i], "-curved"))){
	if(!curved) curved=1;
	randomize=1;
	continue;
      }

      if(!(strcmp(arg_a[i], "-phi") && strcmp(arg_a[i], "-rphi"))){
	rphi=1;
	continue;
      }

      if(!strncmp(arg_a[i], "-curved=", 8)){
	curved=atoi(arg_a[i]+8);
	if(curved<0 || curved>4) curved=1;
	randomize=1;
	continue;
      }

      if(!strncmp(arg_a[i], "-cutth=", 7)){
	test=atof(arg_a[i]+7);
	if(test>=0 && test<=90) cutth=test;
	curved=4;
	randomize=1;
	continue;
      }

      if(!strncmp(arg_a[i], "-EARTHR=", 8)){
	test=atof(arg_a[i]+8);
	if(test>0) EARTHR=test;
	if(!curved) curved=1;
	randomize=1;
	continue;
      }

      if(!strncmp(arg_a[i], "-LENGTH=", 8)){
	test=atof(arg_a[i]+8);
	if(test>0) LENGTH=test;
	randomize=1;
	continue;
      }

      if(!strncmp(arg_a[i], "-RADIUS=", 8)){
	test=atof(arg_a[i]+8);
	if(test>0) RADIUS=test;
	randomize=1;
	continue;
      }

      if(!strncmp(arg_a[i], "-DEPTH=", 7)){
	test=atof(arg_a[i]+7);
	DEPTH=test;
	randomize=1;
	continue;
      }

      if(!strncmp(arg_a[i], "-HEIGHT=", 8)){
	test=atof(arg_a[i]+8);
	if(test>0) HEIGHT=test;
	randomize=1;
	continue;
      }

      if(!strncmp(arg_a[i], "-DCORR=", 7)){
	test=atof(arg_a[i]+7);
	DCORR=test;
	continue;
      }

      if(!strncmp(arg_a[i], "-trigw=", 7)){
	test=atof(arg_a[i]+7);
	trigw=test;
	continue;
      }

      if(!(strcmp(arg_a[i], "-rmpri"))){
	rmpri=1;
	continue;
      }

      if(!(strcmp(arg_a[i], "-rmusr"))){
	rmusr=1;
	continue;
      }

      if(!(strcmp(arg_a[i], "-oms"))){
	oms=1;
	continue;
      }

      if(!(strcmp(arg_a[i], "-ohm"))){
	ohm=1;
	continue;
      }

      if(!(strcmp(arg_a[i], "-msn"))){
	msn=1;
	continue;
      }

      if(!strncmp(arg_a[i], "-cmt=", 5)){
	length=strlen(arg_a[i]+5);
	if(length>=MAX_IN){
	  printf("File name too long\n");
	  return 1;
	}
	if(length<=0){
	  printf("File name too short\n");
	  return 1;
	}
	setcmt=1;
	strcpy(cmtfile, arg_a[i]+5);
	continue;
      }

      if(!strncmp(arg_a[i], "-run=", 5)){
	RUNNUM=atoi(arg_a[i]+5);
	if(RUNNUM<1) RUNNUM=0;
	continue;
      }

      if(!strncmp(arg_a[i], "-FLUXSUM=", 9)){
	FLUXSUM=atof(arg_a[i]+9);
	if(FLUXSUM<=0) FLUXSUM=0;
	continue;
      }

      if(!strncmp(arg_a[i], "-SHOWERS=", 9)){
	SHOWERS=atof(arg_a[i]+9);
	if(SHOWERS<=0) SHOWERS=0;
	continue;
      }

      if(!(strcmp(arg_a[i], "-rr"))){
	remr=1;
	continue;
      }

      if(!strncmp(arg_a[i], "-rr=", 4)){
	REMR=atof(arg_a[i]+4);
	remr=remrflag=1;
	continue;
      }

      if(!(strcmp(arg_a[i], "-test"))){
	if(i+3>=arg_c){
	  printf("Usage: %s -test [num] [theta] [phi]\n", arg_a[0]);
	  return 1;
	}
	test_num=atol(arg_a[i+1]);
	test_cos=atof(arg_a[i+2]);
	test_phi=atof(arg_a[i+3]);
	if(test_num<=0){
	  printf("Number of points must be greater than zero\n");
	  return 1;
	}
	if(test_cos>=90){
	  printf("Unallowed choice of theta; theta must be smaller than 90 degrees\n");
	  return 1;
	}
	if(test_cos<0){
	  printf("Unallowed choice of theta; theta must be greater than zero\n");
	  return 1;
	}
	test_cos*=PI/180;
	test_phi*=PI/180;
	test_cos=cos(test_cos);
	for(i=0; i<test_num; i++){
	  delta(&dx, &dy, test_cos); dx/=test_cos;
	  dx+=DEPTH*sqrt(1/(test_cos*test_cos)-1);
	  test_x=cos(test_phi)*dx-sin(test_phi)*dy;
	  test_y=sin(test_phi)*dx+cos(test_phi)*dy;
	  printf("%g %g\n", test_x, test_y);
	}
	return 0;
      }

      if(!strncmp(arg_a[i], "-cutfe=", 7)){
	cutfe=atof(arg_a[i]+7);
	if(cutfe<=0) cutfe=0;
	continue;
      }

      if(!(strcmp(arg_a[i], "-corr"))){
	rdmcF2k=1;
	continue;
      }

      if(!strncmp(arg_a[i], "-etrans=", 8)){
	etrans=atof(arg_a[i]+8);
	if(etrans<0){
	  printf("Energy of inserted muon must be greater than zero\n");
	  return 1;
	}
	insert=1;
	continue;
      }
      if(!strncmp(arg_a[i], "-dtrans=", 8)){
	dtrans=atof(arg_a[i]+8);
	if(dtrans<0) printf("Separation of inserted muon will be spread randomly between %g and %g\n", -dtrans, gmaxdtrans);
	else if(dtrans==0) printf("Distance of inserted muon must be set. Will use default.\n");
	else printf("Distance of inserted muon will be spread randomly between %g and %g\n", 0.9*dtrans, 1.1*dtrans);
	continue;
      }

      if(!(strcmp(arg_a[i], "-q") && strcmp(arg_a[i], "-quit"))) return 0;
      if(!(strcmp(arg_a[i], "-h") && strcmp(arg_a[i], "-help") && strcmp(arg_a[i], "--help"))){
	printf("\
This is the Utility for CORSIKA output Randomization (ucr), v%s\n\
Available options are: -h  = -help  print this message and quit\n\
                       -q  = -quit  quit immediately\n\
                       -out=[NAME]  change output file to NAME\n\
                       -tr=[NUM]  output only events with # of tracks >= NUM\n\
                       -over=[..] oversample by this number, default is 1\n\
                       -r = -randomize randomize shower core xy locations\n\
                       -c = -curved same for curved Earth's surface\n\
                       -curved=[1-4] different curved surface treatment:\n\
                         1: only downgoing primaries (th from 0 to 90 deg.)\n\
                         2: only upgoing primaries (th from 90 to 180 deg.)\n\
                         3: decide at random, either goes up or down\n\
                         4: oversample x2 each event with th > cutth\n\
                       -phi = -rphi    also randomize azimuth angle\n\
                       -cutth=[degrees]    value used for -curved=4\n\
                       -EARTHR=[radius]    of the Earth            [m]\n\
                       -LENGTH=[length]    of the detector         [m]\n\
                       -RADIUS=[radius]    of the detector         [m]\n\
                       -DEPTH=[depth]      of the detector center  [m]\n\
                       -HEIGHT=[altitude]  of the ice surface      [m]\n\
                       -DCORR=[correction] depth correction (35 m) [m]\n\
                       -trigw=[time in ms] assign event times and combine\n\
                            events that are within trigw ms of each other\n\
                       -rmpri  remove primaries\n\
                       -rmusr  remove user blocks\n\
                       -oms    output only muons\n\
                       -ohm    leave only muon with highest energy per event\n\
                       -msn    leave only muon and neutrinos, delete others\n\
                       -cmt=[file]    append comments contained in the file\n\
                       -run=[number]  set the run number\n\
                       -FLUXSUM=[CORSIKA's value]   per meter2 second sr\n\
                       -SHOWERS=[number] of showers generated by CORSIKA\n\
                       -rr     remove possible previous xy randomization\n\
                       -rr=[DEPTH]   and set the previous value of DEPTH\n\
                       -test [num] [theta] [phi]   test xy randomization\n\
                       -cutfe=[GeV]  angle-dependent cutoff energy for muons\n\
                       -corr   enforce f2k compliance\n\
                       -etrans=[GeV]  insert a muon with given energy and\n\
                       -dtrans=[m]  separation. dtrans<0 sets minimum sep.\n\
", VERSION); return 0; }
      printf("Unknown option: `%s', aborting execution\n", arg_a[i]);
      return 1;
    }
    else flag=1;
  }
  if(flag){
    if(TRNUM==1) printf("All events will output\n");
    else printf("Only events with the number of tracks > %i will output\n", TRNUM);
    if(OVERSAMPLE!=1) printf("Every event will be oversampled %i times\n", OVERSAMPLE);
    if(RUNNUM) printf("Run number will be set to %i\n", RUNNUM);
    if(randomize||rphi||insert){
      srand(RUNNUM);
      printf("Random number seed is initialized with the run number %i\n", RUNNUM);
    }
    if(rphi) printf("Azimuth angle of showers will be randomized\nWarning: magnetic field influence will be destroyed\n");
    if(randomize){
      printf("XY position of shower core location will be randomized\n");
      if(curved){
	printf("Earth's surface will be considered curved (with radius %g km)\n", EARTHR/1.e3);
	if(curved!=4) cutth=0;
	switch(curved){
	case 1: printf("Theta of primaries will range from 0 to 90 degrees\n"); break;
	case 2: printf("Theta of primaries will range from 90 to 180 degrees\n"); break;
	case 3: printf("Theta of primaries will range from 0 to 180 degrees,\n\
events placed in each hemisphere at random\n"); break;
	case 4: printf("Theta of primaries will range from 0 to %g degrees\n\
events with primaries from %g to 90 degrees will be 2x oversampled\n", 180-cutth, cutth);
	cutth=sin(cutth*PI/180)*(1-DEPTH/EARTHR); cutth=sqrt(1-cutth*cutth); break;
	default: printf("Unknown value (%i) of curved.\n", curved); break;
	}
      }
      printf("\tLENGTH = %g\n\tRADIUS = %g\n\tDEPTH  = %g\n\tHEIGHT = %g\n", LENGTH, RADIUS, DEPTH, HEIGHT);
      AREASUM=PI*PI*RADIUS*(RADIUS+LENGTH);
      printf(" AREASUM = %g of meter2 sr\n", AREASUM);
    }
    if(insert){
      if(dtrans<0){
	printf("Inserting a muon with energy of %g GeV and minimum separation of %g m\n\
(minimum transverse momentum of %g GeV/c, assuming %g m int height).\n", etrans, -dtrans, -etrans*dtrans/transheight, transheight);
      }
      else{
	if(dtrans==0){
	  printf("Setting dtrans to default of 400 m\n");
	  dtrans=400.;
	}
	printf("Inserting a muon with energy of %g GeV and random separation between +/- 10 percent of %g m\n\
(minimum transverse momentum of %g GeV/c, assuming %g m int height).\n", etrans, dtrans, 0.9*etrans*dtrans/transheight, transheight);
      }
    }

    if(FLUXSUM!=0) printf(" FLUXSUM = %g per meter2 second sr\n", FLUXSUM);
    LIFETIME=AREASUM*FLUXSUM;
    SHOWERS*=OVERSAMPLE;
    if(LIFETIME!=0 && SHOWERS!=0){
      LIFETIME=SHOWERS/LIFETIME;
      if(curved==3) LIFETIME/=2;
      printf(" %g showers correspond to %g seconds of lifetime\n", SHOWERS, LIFETIME);
    }
    else LIFETIME=0;
    if(trigw>=0){
      if(LIFETIME!=0){
	avtime=LIFETIME;
	if(trigw==0) printf("Assigning event times, events will not be combined\n");
	else{
	  rmusr=1;
	  printf("Assigning event times, event within %g ms will be combined\n", trigw);
	}
      }
      else{
	trigw=-1;
	printf("Cannot compute event times, need valid LIFETIME\n");
      }
    }
    if(DCORR){
      if(randomize) printf("Coordinates will be depth-corrected by %g m\n", DCORR);
      else DCORR=0;
    }
    if(rmpri) printf("Primary tracks will be removed\n");
    if(rmusr) printf("User blocks will be removed\n");
    if(ohm) oms=1;
    if(msn) oms=0, ohm=0;
    if(oms) printf("Only muons will output\n");
    if(ohm) printf("Only highest energy muons will output\n");
    if(msn) printf("Only muons and neutrinos will output\n");
    if(remr){
      printf("Possible previous xy randomization will be removed\n");
      if(!remrflag) printf("Warning: this option may give incorrect results when re-randomizing.\n\
Please specify the previous value of DEPTH explicitly by -rr=[DEPTH]\n");
      else printf("Previous value of DEPTH is set to %g\n", REMR);
    }
    if(cutfe!=0){
      llow=log(1+(cutfe*blow/alow))/blow;
      printf("Skew angle cutoff at %g GeV for the level of %g mwe is enabled\n", cutfe, llow);
    }
    if(setcmt) printf("Comments file `%s' will be appended\n", cmtfile);
    if(rdmcF2k) printf("Enforcing F2k compliance\n");
    printf("Output file name is `%s'\n", outfile);
    printf("\nProcessing files:\n");
    if((fo=rdmc_mcopen(outfile, "w", F2000_F))==NULL){
      printf("cannot open output file `%s', exiting\n", outfile);
      return 1;
    }

    if(stat(outfile, &outstat)) exit(1);
    rdmc_write_parameters(fo, arg_c, arg_a, VERSION);
    if(FLUXSUM){
      sprintf(s, " ! FLUXSUM is %g per meter2 second sr\n", FLUXSUM);
      rdmc_concat_comment(&(fo->creator), s, F2000_F);
    }

    if(AREASUM){
      sprintf(s, " ! AREASUM is %g of meter2 sr\n", AREASUM);
      rdmc_concat_comment(&(fo->creator), s, F2000_F);
    }

    if(LIFETIME){
      sprintf(s, " ! LIFETIME per %g showers is %g seconds\n", SHOWERS, LIFETIME);
      rdmc_concat_comment(&(fo->creator), s, F2000_F);
    }

    if(setcmt){
      if((fc=fopen(cmtfile, "r"))==NULL){
	printf("cannot open comments file `%s', exiting\n", cmtfile);
	return 1;
      }
      rdmc_concat_comment(&(fo->creator), " ! Card file follows:\n", F2000_F);
      while(!my_input(s, fc)) rdmc_concat_comment(&(fo->creator), s, F2000_F);
      fclose(fc);
    }

    if(avtime>0){
      jmax=OVERSAMPLE;
      OVERSAMPLE=1;
      if(curved==4) printf("Warning: using `-curved=4' together with -trigw=%g\n\
  may lead to undefined results\n", trigw);
    }
    else jmax=1;

    for(j=0; j<jmax; j++) for(i=1; i<arg_c; i++) if(*arg_a[i]!='-'){

      if(avtime>0){
 	printf("analyzing file `%s'\n", arg_a[i]);
 	if((fi=rdmc_mcopen(arg_a[i], "r", F2000_F))==NULL){
	  printf("cannot open input file `%s', exiting\n", arg_a[i]);
	  return 1;
	}
	rdmc_rarr(fi, &ar);
	evts=0;
	while(1){
	  if(rdmc_revt(fi, &ar, &ev)) break;
	  evts++;
	}
	rdmc_mcclose(fi);
	printf("%li events in this file\n", evts);
	if(evts>0) avtime=LIFETIME/(jmax*evts);
	else{
	  avtime=0;
	  printf("File %s contains no usable events\n", arg_a[i]);
	}
      }

      printf("file `%s' is being fed in\n", arg_a[i]);
      errf=process(arg_a[i], fo);
      if(errf) printf("error in function `process'\n");
      error+=errf;
    }
    rdmc_mcclose(fo);
  }
  else printf("You must specify at least one input file: `%s [infile]'\n", arg_a[0]);

  return error;
}