static const float xx=1.e-5;
static const float zoff=1948.07;

static const float OVR=5;
static const float R=OVR*0.178;

static const float ppm=2107.84/(OVR*OVR);
static const float crf[4]={ ppm*4.889*0.894, // em shower
			    ppm*4.076*0.860, // hdr m/GeV
			    ppm*1.1720,      // bare
			    ppm*0.0324 };    // muon

static const int wnum=40;
struct wv{
  float ocm;                  // 1 / speed of light in medium
  float coschr, sinchr;       // cos and sin of the cherenkov angle
  float *abs, *sca;           // arrays of absorption and scattering values
} wvd[wnum];

struct ice {
  wv *w;
  float g;                    // <cos(scattering angle)>
  float ocv;                  // 1 / speed of light in vacuum
  int size;                   // size of kurt table
  float dh, hdh, rdh, hmin;   // step, step/2, 1/step, and min depth

  ice(){   // initialize ice parameters
    ifstream inFile;

    bool flag=true;
    inFile.open("wv.dat", ifstream::in);
    vector<float> wx, wy;

    if(!inFile.fail()){
      float xa, ya, xo=0, yo=0;
      int num=0;
      while(inFile>>xa>>ya){
	if(( xa<0 || 1<xa ) || num==0 && xa!=0 || num>0 && ( xa<=xo || ya<=yo )){ flag=false; break; }
	wx.push_back(xa); wy.push_back(ya);
	xo=xa; yo=ya; num++;
      }
      if(xo!=1 || wx.size()<2) flag=false;
      inFile.close();
      if(flag){ cerr<<"Loaded "<<wx.size()<<" wavelenth points"<<endl; }
      else{ cerr<<"File wv.dat did not contain valid data"<<endl; exit(1); }
    }
    else{ cerr<<"Could not open file wv.dat"<<endl; exit(1); }

    flag=true;
    float A, B, D, E, a, k;
    float Ae, Be, De, Ee, ae, ke;
    vector<float> dp, be, ba, td;

    inFile.open("icemodel.par", ifstream::in);

    if(flag=!inFile.fail()){
      if(flag) flag=(inFile >> a >> ae);
      if(flag) flag=(inFile >> k >> ke);
      if(flag) flag=(inFile >> A >> Ae);
      if(flag) flag=(inFile >> B >> Be);
      if(flag) flag=(inFile >> D >> De);
      if(flag) flag=(inFile >> E >> Ee);
      if(!flag) cerr << "File icemodel.par found, but is corrupt" << endl;
      inFile.close();
      if(!flag) exit(1);
    }
    else{ cerr << "File icemodel.par was not found" << endl; exit(1); }

    inFile.open("icemodel.dat", ifstream::in);
    if(!inFile.fail()){
      size=0;
      float dpa, bea, baa, tda;
      while(inFile >> dpa >> bea >> baa >> tda){
	dp.push_back(dpa);
	be.push_back(bea);
	ba.push_back(baa);
	td.push_back(tda);
	size++;
      }
      inFile.close();
      if(size<2){ cerr << "File icemodel.dat found, but is corrupt" << endl; exit(1); }
    }
    else{ cerr << "File icemodel.dat was not found" << endl; exit(1); }

    dh=dp[1]-dp[0];
    if(dh<=0){ cerr << "Ice table does not use increasing depth spacing" << endl; exit(1); }

    for(int i=0; i<size; i++) if(i>0) if(fabs(dp[i]-dp[i-1]-dh)>dh*xx){
      cerr << "Ice table does not use uniform depth spacing" << endl; exit(1);
    }
    cerr<<"Loaded "<<size<<" ice layers"<<endl;

    hdh=dh/2; rdh=1/dh; hmin=zoff-dp[size-1];
    g=0.8;

    int nfla;
    float wfla=405, fmin;
    for(int n=0; n<wnum; n++){
      float p=(n+0.5)/wnum;
      int m=0; while(wx[++m]<p);
      float wva=(wy[m-1]*(wx[m]-p)+wy[m]*(p-wx[m-1]))/(wx[m]-wx[m-1]);

      float faux=abs(wva-wfla);
      if(n==0 || faux<fmin){ nfla=n; fmin=faux; }

      float wv0=400;
      float l_a=pow(wva/wv0, -a);
      float l_k=pow(wva, -k);
      float ABl=A*exp(-B/wva);

      w=&wvd[n];
      w->abs = new float[size];
      w->sca = new float[size];

      for(int i=0; i<size; i++){
	int j=size-1-i;
	w->sca[i]=be[j]*l_a/(1-g); w->abs[i]=(D*ba[j]+E)*l_k+ABl*(1+0.01*td[j]);
	if(w->sca[i]<=0 || w->abs[i]<=0){ cerr << "Invalid value of ice parameter, cannot proceed" << endl; exit(1); }
      }

      float wv=wva*1.e-3;
      float wv2=wv*wv;
      float wv3=wv*wv2;
      float wv4=wv*wv3;
      float np=1.55749-1.57988*wv+3.99993*wv2-4.68271*wv3+2.09354*wv4;
      float ng=np*(1+0.227106-0.954648*wv+1.42568*wv2-0.711832*wv3);
      float c=0.299792458; ocv=1/c; w->ocm=ng/c;
      w->coschr=1/np; w->sinchr=sqrt(1-w->coschr*w->coschr);
    }
    w=&wvd[nfla];
  }

  void destroy(){   // empties booked arrays and tables
    for(int n=0; n<wnum; n++){
      delete wvd[n].abs;
      delete wvd[n].sca;
    }
  }

  ~ice(){
    destroy();
  }

  void next(float rnd){
    w=&wvd[lroundf(wnum*rnd-0.5)];
  }

  // get distance for overburden x
  float aX(float x, float z, float nz){
    int i=lroundf((z-hmin)*rdh), j;
    if(i<0) i=0; else if(i>=size) i=size-1;

    float h=hmin+i*dh;
    float ax, ai, ah;

    if(nz>0){
      ah=h+hdh;
      ai=(nz*x-(ah-z)*w->abs[i])*rdh;
      for(j=i; j<size-1 && ai>0; ah+=dh) ai-=w->abs[++j];
    }
    else{
      ah=h-hdh;
      ai=(nz*x-(ah-z)*w->abs[i])*rdh;
      for(j=i; j>0 && ai<0; ah-=dh) ai+=w->abs[--j];
    }

    if(i==j) ax=x/w->abs[i];
    else ax=(ai*dh/w->abs[j]+ah-z)/nz;

    return ax;
  }

  // get distance for overburden x
  float sX(float x, float z, float nz){
    int i=lroundf((z-hmin)*rdh), j;
    if(i<0) i=0; else if(i>=size) i=size-1;

    float h=hmin+i*dh;
    float sx, ai, ah;

    if(nz>0){
      ah=h+hdh;
      ai=(nz*x-(ah-z)*w->sca[i])*rdh;
      for(j=i; j<size-1 && ai>0; ah+=dh) ai-=w->sca[++j];
    }
    else{
      ah=h-hdh;
      ai=(nz*x-(ah-z)*w->sca[i])*rdh;
      for(j=i; j>0 && ai<0; ah-=dh) ai+=w->sca[--j];
    }

    if(i==j) sx=x/w->sca[i];
    else sx=(ai*dh/w->sca[j]+ah-z)/nz;

    return sx;
  }

  // get overburden for distance a
  float xA(float a, float z, float nz){
    int i=lroundf((z-hmin)*rdh), j;
    if(i<0) i=0; else if(i>=size) i=size-1;

    float h=hmin+i*dh;
    float y=z+nz*a;

    j=lroundf((y-hmin)*rdh);
    if(j<0) j=0; else if(j>=size) j=size-1;

    float xa;

    if(i==j) xa=a*w->abs[i]; 
    else{
      float g=hmin+j*dh;
      if(nz>0){
	xa=(y-(g-hdh))*w->abs[j];
	while(--j>i) xa+=dh*w->abs[j];
	xa+=(h+hdh-z)*w->abs[i];
      }
      else{
	xa=(y-(g+hdh))*w->abs[j];
	while(++j<i) xa-=dh*w->abs[j];
	xa+=(h-hdh-z)*w->abs[i];
      }
      xa/=nz;
    }

    return xa;
  }
} m;