// #define STRETCH

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, rdh, hmin, hmax;  // step, 1/step, min and max depth
  double *abs, *sca;           // arrays of absorption and scattering values
  double *rabs, *rsca;         // inverse of these arrays
  double *abs_sum, *sca_sum;   // arrays of absorption and scattering depth-integrated values
  bool kurt;                   // whether kurt table is used
  bool accl;                   // accelerate lookup for large ice tables

  bool first;
  double A, B, D, E, a, k;
  double Ae, Be, De, Ee, ae, ke;
  vector<double> dp, be, ba, td;

#ifdef STRETCH
  double stretch(double x, double mean, double factor){
    return mean+factor*(x-mean);
  }
#endif

  void init(){  // initialize ice parameters
    first=true;
    kurt=false;

    bool flag=true;
    ifstream inFile;
    inFile.open((dir+"icemodel.par").c_str(), 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 return;

    inFile.open((dir+"icemodel.dat").c_str(), ifstream::in);
    if(!inFile.fail()){
      size=0;
      double dpa, bea, baa, tda;
      while(inFile >> dpa >> bea >> baa >> tda){
#ifdef STRETCH
	if(dpa<2090){
	  bea=stretch(bea, 0.047, 2);
	  baa=stretch(baa, 0.041, 2);
	}
	else{
	  bea=stretch(bea, 0.029, 2);
	  baa=stretch(baa, 0.017, 2);
	}
#endif
	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;

    dh=dp[1]-dp[0]; rdh=1/dh;
    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; }
    }

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

    kurt=true;
  }

  void init(double wv){
    this->wv=wv;

    if(kurt){
      if(wv<0.2) wv=0.2;
      else if(wv>0.8) wv=0.8;
      double wva=wv*1.e3;

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

      if(first){
	abs = new double[size];
	sca = new double[size];
	rabs = new double[size];
	rsca = new double[size];
      }

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

      if(first){
	sca_a=size/sca_a;
	abs_a=size/abs_a;

	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;
	rabs[i]=1/abs[i];
	rsca[i]=1/sca[i];
	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);

      if(first) 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;
    }
    else{
      if(first) cerr << "Using bulk ice model at l=" << wv*1.e3 << "nm: np=" << np << " ng=" << ng << " sca=" << le << " abs=" << la << endl;
    }
    first=false;
  }

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

  ~ice(){
    destroy();
  }

  // constructor
  ice(){
    accl=false;

    c=0.299792458; g=0.8;
    la=105.; le=21.8;
    ls=(1-g)*le;
    wv=0.380;

    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();
    init(wv);
  }

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

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

  // 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)*rdh); 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)*rdh;
    int j=i;
    if(nz<0) getlo(j, abs_sum, dif); else gethi(j, abs_sum, dif);
    if(i==j) tot=x*rabs[i];
    else tot=((hmin+(j+1)*dh)-(abs_sum[j]-dif)*dh*rabs[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)*rdh); 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)*rdh); 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)*rdh); 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)*rdh;
    int j=i;
    if(nz<0) getlo(j, sca_sum, dif); else gethi(j, sca_sum, dif);
    if(i==j) tot=x*rsca[i];
    else tot=((hmin+(j+1)*dh)-(sca_sum[j]-dif)*dh*rsca[j]-z)/nz;
    return tot;
  }
} m;