struct ice {
  double la;    // absorption length
  double ls;    // scattering length
  double le;    // effective scattering length
  double g;     // <cos(scattering angle)>

  double wv;    // light wavelength
  double np;    // phase refractive index
  double ng;    // group refractive index

  double c;     // speed of light in vacuum
  double cm;    // speed of light in medium
  double na;    // average propagation refractive index
  double coschr, sinchr;       // cos and sin of the cherenkov angle

  int size;                    // size of kurt table
  double dh, hmin, hmax;       // step, min and max depth
  double *abs, *sca;           // arrays of absorption and scattering values
  double *abs_sum, *sca_sum;   // arrays of absorption and scattering depth-integrated values
  bool kurt;                   // whether kurt table is used

  void init(){  // initialize ice parameters
    if(wv<0.3) wv=0.3;
    else if(wv>0.6) wv=0.6;
    double wva=wv*1.e3;

    kurt=false;

    double A, B, D, E, a, k;
    double Ae, Be, De, Ee, ae, ke;

    bool flag=true;
    ifstream inFile;
    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, falling back to bulk ice" << endl;
      inFile.close();
      if(!flag) return;
    }
    else{
      cerr << "Using bulk ice model at l=" << wva << "nm: np=" << np << " ng=" << ng
	   << " sca=" << le << " abs=" << la << endl;
      return;
    }

    vector<double> dp, be, ba, td;
    inFile.open("icemodel.dat", ifstream::in);
    if(!inFile.fail()){
      size=0;
      double 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++;
      }
      if(size<2){ cerr << "File icemodel.dat found, but is corrupt, falling back to bulk ice" << endl; return; }
      inFile.close();
    }
    else return;

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

    abs = new double[size];
    sca = new double[size];

    double abs_a=0, sca_a=0;
    dh=dp[1]-dp[0];
    if(dh<=0){ cerr << "Ice table does not use increasing depth spacing, falling back to bulk ice" << endl; return; }

    for(int i=0; i<size; i++){
      if(i>0) if(fabs(dp[i]-dp[i-1]-dh)>dh*1.e-10){ cerr << "Ice table does not use uniform depth spacing, falling back to bulk ice" << endl; return; }
      sca[i]=be[i]*l_a;
      abs[i]=(D*ba[i]+E)*l_k+ABl*(1+0.012*td[i]);
      if(sca[i]<=0 || abs[i]<=0){ cerr << "Invalid value of ice parameter, falling back to bulk ice" << endl; return; }
      sca_a+=sca[i]; abs_a+=abs[i];
    }

    sca_a=size/sca_a;
    abs_a=size/abs_a;

    hmin=dp[0]-dh/2;
    hmax=dp[size-1]+dh/2;

    abs_sum = new double[size];
    sca_sum = new double[size];
    double abs_aux=0, sca_aux=0;

    for(int i=0; i<size; i++){
      sca[i]/=1-g;
      abs_aux+=abs[i];
      sca_aux+=sca[i];
      abs_sum[i]=abs_aux;
      sca_sum[i]=sca_aux;
    }

    double wv2=wv*wv;
    double wv3=wv*wv2;
    double wv4=wv*wv3;
    np=1.55749-1.57988*wv+3.99993*wv2-4.68271*wv3+2.09354*wv4;
    ng=np*(1+0.227106-0.954648*wv+1.42568*wv2-0.711832*wv3);
    na=(np*ng-1)/sqrt(np*np-1); cm=c/ng;
    coschr=1/np; sinchr=sqrt(1-coschr*coschr);

    cerr << "Using 6-parameter ice model at l=" << wva << "nm: np=" << np << " ng=" << ng
	 << " average sca=" << sca_a << " abs=" << abs_a << " depths=(" << hmin << ";" << hmax << ")" << endl;
    kurt=true;
  }

  void destroy(){  // empties booked arrays and tables
    if(kurt){
      delete abs;
      delete sca;
      delete abs_sum;
      delete sca_sum;
    }
  }

  ~ice(){
    destroy();
  }

  // constructor
  ice(){
    c=0.299792458; g=0.8;
    la=105.; le=21.8;
    ls=(1-g)*le;
    wv=0.405;

    np=1.31950; ng=1.35634;
    na=(np*ng-1)/sqrt(np*np-1); cm=c/ng;
    coschr=1/np; sinchr=sqrt(1-coschr*coschr);

    init();
  }

  void gethi(int & j, double a[], double max){
    int j1=j, j2=size-1;
    double f1=a[j1], f2=a[j2];
    while(j2-j1>1){
      int jc=int((j1+j2)/2);
      double fa=a[jc];
      if(fa>max) j2=jc; else j1=jc;
    }
    j=j1;
  }

  void getlo(int & j, double a[], double min){
    int j1=0, j2=j;
    double f1=a[j1], f2=a[j2];
    while(j2-j1>1){
      int jc=int((j1+j2)/2);
      double fa=a[jc];
      if(fa>min) j2=jc; else j1=jc;
    }
    j=j2;
  }

  // get distance for overburden x
  double aX(double x, double z, double nz){
    if(!kurt) return x*la;
    z=-z; nz=-nz;
    int i=(int) floor((z-hmin)/dh); if(i<0) i=0; else if(i>=size) i=size-1;
    double tot=0;

    double ai=((hmin+(i+1)*dh)-z)*abs[i];
    double dif=abs_sum[i]+(x*nz-ai)/dh;
    int j=i;
    if(nz<0){ getlo(j, abs_sum, dif); for(; j>0; j--) if(abs_sum[j-1]<dif) break; }
    else{ gethi(j, abs_sum, dif); for(; j<size-1; j++) if(abs_sum[j]>dif) break; }
    if(i==j) tot=x/abs[i];
    else tot=((hmin+(j+1)*dh)-(abs_sum[j]-dif)*dh/abs[j]-z)/nz;
    return tot;
  }

  // get overburden for distance a
  double xA(double a, double z, double nz){
    if(!kurt) return a/la;
    z=-z; nz=-nz;
    int i=(int) floor((z-hmin)/dh); if(i<0) i=0; else if(i>=size) i=size-1;
    double tot=0;

    double y=z+nz*a;
    int j=(int) floor((y-hmin)/dh); if(j<0) j=0; else if(j>=size) j=size-1;
    if(i==j) tot=a*abs[i];
    else tot=((abs_sum[j]-abs_sum[i])*dh-((hmin+(j+1)*dh)-y)*abs[j]+((hmin+(i+1)*dh)-z)*abs[i])/nz;
    return tot;
  }

  // get distance for overburden x
  double sX(double x, double z, double nz){
    if(!kurt) return x*ls;
    z=-z; nz=-nz;
    int i=(int) floor((z-hmin)/dh); if(i<0) i=0; else if(i>=size) i=size-1;
    double tot=0;

    double ai=((hmin+(i+1)*dh)-z)*sca[i];
    double dif=sca_sum[i]+(x*nz-ai)/dh;
    int j=i;
    if(nz<0){ for(; j>0; j--) if(sca_sum[j-1]<dif) break; }
    else{ for(; j<size-1; j++) if(sca_sum[j]>dif) break; }
    if(i==j) tot=x/sca[i];
    else tot=((hmin+(j+1)*dh)-(sca_sum[j]-dif)*dh/sca[j]-z)/nz;
    return tot;
  }
} m;