dats d;

static const float zoff=1948.07;
static const int wnum=40;
ices * iced[wnum];

struct ini{
  ini(){

    { // initialize random numbers
      int size;

      ifstream inFile;
      inFile.open("rnd.txt", ifstream::in);
      vector<unsigned int> rx;

      if(!inFile.fail()){
	unsigned int a;
	unsigned long long b,c;
	while(inFile>>a>>b>>c) rx.push_back(a);
	if(rx.size()<1){ cerr<<"File rnd.txt did not contain valid data"<<endl; exit(1); }
      }
      else{ cerr<<"Could not open file rnd.txt"<<endl; exit(1); }
      size=rx.size();
      cerr<<"Loaded "<<size<<" random multipliers"<<endl;

      d.rm = new unsigned int[size];
      d.rs = new unsigned long long[size];

      state s;
      s.rsize=size;
      s.ml=2083801278;
      s.ax=362436069;
      s.dx=1234567;

      for(int i=0; i<size; i++){
	d.s=s;
	d.rm[i]=rx[i];
	d.rs[i]=s.da;
      }
    }

    { // initialize ice parameters
      ices *w;
      int size;                   // size of kurt table
      float g;                    // <cos(scattering angle)>
      float dh, hdh, rdh, hmin;   // step, 1/step, min and max depth

      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(fabsf(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=fabsf(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 = (ices *) new float[2*size+sizeof(ices)/sizeof(float)];
	iced[n] = w;

	d.g=g;
	d.size=size;
	d.dh=dh;
	d.hdh=hdh;
	d.rdh=rdh;
	d.hmin=hmin;

	for(int i=0; i<size; i++){
	  int j=size-1-i;
	  w->z[i].sca=be[j]*l_a/(1-g); w->z[i].abs=(D*ba[j]+E)*l_k+ABl*(1+0.01*td[j]);
	  if(w->z[i].sca<=0 || w->z[i].abs<=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; d.ocv=1/c; w->ocm=ng/c;
	w->coschr=1/np; w->sinchr=sqrt(1-w->coschr*w->coschr);
      }
      d.w=iced[nfla];
    }

    { // initialize geometry
      ifstream inFile;
      inFile.open("geo-f2k", ifstream::in);
      vector<DOM> oms;

      if(!inFile.fail()){
	DOM dom;
	string mbid;
	int om, str;
	unsigned long long omid;
	while(inFile>>mbid>>hex>>omid>>dec>>dom.r[0]>>dom.r[1]>>dom.r[2]>>str>>om) if(str>0 && om>=1 && om<=60){
	  dom.r[2]+=zoff;
	  dom.n.dom=om;
	  dom.n.str=str;
	  oms.push_back(dom);
	}
	inFile.close();
      }

      d.gsize = oms.size();
      d.oms = new DOM[d.gsize];
      for(int n=0; n<d.gsize; n++) d.oms[n]=oms[n];

#ifdef ACCL2
      d.cn[0]=cx;
      d.cn[1]=cy;
      d.cn[2]=cz;
      float cl[3], ch[3], crst[3];

      vector<int> cells[cx][cy][cz];
      for(int i=0; i<cx; i++) for(int j=0; j<cy; j++) for(int k=0; k<cz; k++)
	cells[i][j][k].erase(cells[i][j][k].begin(), cells[i][j][k].end());

      for(int n=0; n<oms.size(); n++) for(int m=0; m<3; m++){
	if(n==0 || oms[n].r[m]<cl[m]) cl[m]=oms[n].r[m];
	if(n==0 || oms[n].r[m]>ch[m]) ch[m]=oms[n].r[m];
      }
      for(int m=0; m<3; m++){
	float diff=ch[m]-cl[m];
	if(diff<R){ cerr<<"Error in cell initialization"<<endl; exit(1); }
	crst[m]=(d.cn[m]-1)/diff;
      }

      for(int n=0; n<oms.size(); n++){
	int i[3];
	for(int m=0; m<3; m++){
	  i[m]=lroundf((oms[n].r[m]-cl[m])*crst[m]);
	  if(i[m]<0 && i[m]>=d.cn[m]){ cerr<<"Error in cell initialization"<<endl; exit(1); }
	}
	cells[i[0]][i[1]][i[2]].push_back(n);
      }

      for(int m=0; m<3; m++){ d.cl[m]=cl[m]; d.crst[m]=crst[m]; }
      d.carr = new short[cx*cy*cz+d.gsize];
      int m=0;
      for(int i=0; i<cx; i++) for(int j=0; j<cy; j++) for(int k=0; k<cz; k++){
	d.cidx[i][j][k]=m;
	for(int l=0; l<cells[i][j][k].size(); l++, m++) d.carr[m]=cells[i][j][k][l];
	d.carr[m++]=-1;
      }
#endif

      cerr<<"Loaded "<<oms.size()<<" DOMs"<<endl;
    }

    {
      d.hidx=0;
      d.hnum=100000;
      d.hits = new hit[d.hnum];
    }
  }

  ~ini(){   // empties booked arrays and tables
    delete d.rm;
    delete d.rs;
    for(int n=0; n<wnum; n++) delete iced[n];
    delete d.oms;
#ifdef ACCL2
    delete d.carr;
#endif
    delete d.hits;
  }

  void next(float rnd){
    d.w=iced[lroundf(wnum*rnd-0.5)];
  }
} m;