/*======================================================================
 * numi-flux:
 * ==========
 * [0] Compile with:
 *   gcc numi-flux.c -o numi-flux -lm
 *
 * [1] Set the environment variable NUMI_FLUX_DIR to the directory where
 * the flux text files are located. The flux files can be downloaded
 * from:
 *   http://enrico1.physics.indiana.edu/messier/off-axis/numi-fluxes.tar.gz
 *
 * [2] Run the program with:
 *
 *   ./numi-flux [tune] [x] [y] [z] > my-flux-table.txt
 *
 * where [tune] is a NuMI horn tune. For example, le for the low
 * energy beam tune, me for the medium energy tune, lerev for the le
 * tune with horn current reversed, and merev for the medium energy
 * tune with horn currents reversed. The detector position should be
 * specified in units of km. The system of units has the z-axis along
 * the NuMI decay tunnel with y up and x pointing to beam left.
 * 
 * [3] The flux tables produced contain seven columns:
 *      1 - neutrino energy in GeV at bin center
 *      2 - the nue flux
 *      3 - the numu flux
 *      4 - the nutau flux
 *      5 - the nue-bar flux
 *      6 - the numu-bar flux
 *      7 - the nutau-bar flux
 *    The units for fluxes are #neutrinos/cm^2/bin/1E20POT. The tables
 *    are made in 1000 bins of 50 MeV each ranging from 0 GeV to 50
 *    GeV.
 *
 * [4] The program uses a set of fluxes produced at a fixed z location
 * (1000 km) and various off-axis locations ranging from 0 to 100
 * km. The fluxes at new locations are interpolated linearly and then
 * scaled by 1/r^2. I've spot checked several sites and the
 * interpolated results agree typically within 10% of what I get
 * calculating the sites in detail. As such, these results should be
 * used to "explore" the fluxes. Once you settle on a site, the fluxes
 * should be calcualted in detail.
 *
 * [5] Please send bug reports to messier@indiana.edu
 *======================================================================
 */
#include "stdio.h"
#include "stdlib.h"
#include "math.h"

static float gsenergy[1000];    /* Energy at center of bin (GeV) */
static float gsfluxlo[1000][6]; /* Fluxes at lower angle limit */
static float gsfluxhi[1000][6]; /* Fluxes at upper angle limit */
static float gsflux[1000][6];   /* Interpolated fluxes */

void numiflx_interp(const char* tune, float x, float y, float z) 
{
  /*====================================================================
   * Fill arrays with interpolated fluxes. 
   *
   * Inputs:
   *   tune  - the NuMI target/horn tune (le, lerev, me, merev)
   *   x,y,z - the detector location in units of km
   *
   * Outputs:
   *   Written to the gsflux[][] array
   *
   * The routine uses flux tables pre-calculated at z=1000 km and
   * x=0,100 km in steps of 2 km. Based on the input detector
   * location, the two flux tables which bracket the input location
   * are loaded. The fluxes are then linearly interpolated to account
   * for the difference in the detector's angular location. Finally,
   * the fluxes are scaled by 1/r^2 to account for the different
   * detector distances from Fermilab.
   ====================================================================*/
  int         i;          /* counter */
  int         j;          /* counter */
  int         ilo;        /* x position of lower angle bracket */
  int         ihi;        /* x position of upper angle bracket */
  float       ang;        /* Angle of requested detector location */
  float       anglo;      /* Angle of lower bracket */
  float       anghi;      /* Angle of upper bracker */
  float       z2 =1000.0; /* z location of template fluxes */
  float       t2;         /* Position of requested location projected to z2 */
  float       r;          /* Distance to requested location */
  float       rlo;        /* Distance to location in lower bracket */
  float       rhi;        /* Distance to location in upper bracket */
  char        file1[256]; /* File name for lower bracket flux table */
  char        file2[256]; /* File name for upper bracket flux table */
  char        buf[256];   /* File buffer */
  const char* dir = 0;    /* Directory name for flux tables */
  FILE*       fp;         /* Flux file pointer */
  
  /* Work out geometry of requested location */
  ang = atan2(sqrt(x*x + y*y),sqrt(x*x + y*y + z*z));
  t2  = z2*sin(ang);
  r   = sqrt(x*x + y*y + z*z);

  /* Flux files are calculated every 2 km at z = 1000 km. Find next
     lowest and highest transverse positions which are even number of
     km. These locations bracket the requested location */
  ilo = (int)rint(t2);
  if (ilo%2 != 0) --ilo;
  ihi = ilo+2;

  /* Work out the distances and angles for each of the two flux files */
  anglo = atan2((float)ilo, z2);
  anghi = atan2((float)ihi, z2);
  rlo   = sqrt((float)(ilo*ilo)+z2*z2);
  rhi   = sqrt((float)(ihi*ihi)+z2*z2);

  /* Find the directory for the flux files and make the filenames */
  dir = getenv("NUMI_FLUX_DIR");
  if (dir==0) {
    fprintf(stderr,
	    "numi-flux.c: "
	    "Set NUMI_FLUX_DIR to point to location of the "
	    "NUMI flux files\n");
    exit(1);
  }
  sprintf(file1,"%s/%s-%d-%d.txt",dir, tune, 1000, ilo);
  sprintf(file2,"%s/%s-%d-%d.txt",dir, tune, 1000, ihi);
  
  /* Open the file for the lower angle and read it in */
  fp = fopen(file1,"r");
  if (fp==0) {
    fprintf(stderr,"numi-flux.c: Cannot open file %s for read.\n",
	    file1);
    exit(1);
  }
  for (i=0; i<1000; ++i) {
    fgets(buf,256,fp);
    sscanf(buf,"%f %f %f %f %f %f %f\n",
	   &gsenergy[i],
	   &gsfluxlo[i][0], &gsfluxlo[i][1], &gsfluxlo[i][2],
	   &gsfluxlo[i][3], &gsfluxlo[i][4], &gsfluxlo[i][5]);
    for (j=0; j<6; ++j) gsfluxlo[i][j] *= rlo*rlo;
  }
  close(fp);

  /* Open the file for the upper angle and read it in */
  fp = fopen(file2,"r");
  if (fp==0) {
    fprintf(stderr,"numi-flux.c: Cannot open file %s for read.\n",
	    file2);
    exit(1);
  }
  for (i=0; i<1000; ++i) {
    fgets(buf,256,fp);
    sscanf(buf,"%f %f %f %f %f %f %f\n",
	   &gsenergy[i],
	   &gsfluxhi[i][0], &gsfluxhi[i][1], &gsfluxhi[i][2],
	   &gsfluxhi[i][3], &gsfluxhi[i][4], &gsfluxhi[i][5]);
    for (j=0; j<6; ++j) gsfluxhi[i][j] *= rhi*rhi;
  }
  close(fp);

  /* Make the interpolation in off-axis angle and scale by 1/r^2 */
  for (i=0; i<1000; ++i) {
    for (j=0; j<6; ++j) {
      gsflux[i][j] =
	gsfluxlo[i][j] +
	(ang-anglo)*(gsfluxhi[i][j]-gsfluxlo[i][j])/(anghi-anglo);
      gsflux[i][j] /= (r*r);
    }
  }
  
  /* Dump the table to the screen */
  for (i=0; i<1000; ++i) {
    printf("%7.3f %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e\n",
	   gsenergy[i],
	   gsflux[i][0],
	   gsflux[i][1],
	   gsflux[i][2],
	   gsflux[i][3],
	   gsflux[i][4],
	   gsflux[i][5]);
  }
}

/**********************************************************************/
int main(int argc, char** argv) 
{
  if (argc==1 || strcmp(argv[1],"-h")==0) {
    printf("Usage: numi-flux [tune] [x] [y] [z]\n"
	   "  [tune] = NuMI horn tune (le, me, lerev, merev)\n"
	   "  [x]    = Detector x position (km)\n"
	   "  [y]    = Detector x position (km)\n"
	   "  [z]    = Detector x position (km)\n");
    exit(0);
  }
  numiflx_interp(argv[1],atof(argv[2]),atof(argv[3]),atof(argv[4]));
  return 0;
}

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