const float fcv=FPI/180.f; static unsigned int ovr=OVER; int xr=1; struct ikey{ int str, dom; bool isinice() const{ return str>0 && dom>=1 && dom<=60; } bool operator< (const ikey & rhs) const { return str == rhs.str ? dom < rhs.dom : str < rhs.str; } bool operator!= (const ikey & rhs) const { return str != rhs.str || dom != rhs.dom; } }; struct OM:DOM,ikey{}; vector i3oms; struct name:ikey{ int type; float rde, hv; name(){} name(ikey k, float r, int t, float h):ikey(k){ rde=r; type=t; hv=h; } }; map hvs; map > rdes; bool rdef=false; float rmax=0, rmay=0; dats d; datz z; struct doms{ DOM oms[MAXGEO]; name names[MAXGEO]; map rde; hit * hits; photon * pz; #ifdef TALL pbuf * bf; #endif } q; unsigned char sname(int n){ name & s = q.names[n]; return s.str>78&&s.dom>10?s.str+10:s.str; } static const float zoff=1948.07; unsigned int sv=0; void rs_ini(){ union{ unsigned long long da; struct{ unsigned int ax; unsigned int dx; }; } s; s.ax=362436069; s.dx=1234567; s.ax+=sv; for(int i=0; i v; while(getline(inFile, in)) if(sscanf(in.c_str(), "%f", &aux)==1) v.push_back(aux); if(v.size()>=4){ int xR=lroundf(v[0]); d.zR=1.f/xR, ovr*=xR*xR, d.R=OMR*xR; d.R2=d.R*d.R; d.eff=v[1], d.sf=v[2], d.g=v[3]; d.g2=d.g*d.g; d.gr=(1-d.g)/(1+d.g); xr=xR; cerr<<"Configured: xR="<=7){ const float cv=FPI/180, thx=v[4]; d.azx=cos(cv*thx), d.azy=sin(cv*thx); d.k1=exp(v[5]); d.k2=exp(v[6]); d.kz=1/(d.k1*d.k2); d.k=isfinite(d.kz) && ( d.k1!=1 || d.k2!=1 )?1:0; if(d.k>0) cerr<<"Ice anisotropy is k("< rx; ifstream inFile((dir+"rnd.txt").c_str(), ifstream::in); if(!inFile.fail()){ string in; while(getline(inFile, in)){ stringstream str(in); unsigned int a; if(str>>a) rx.push_back(a); } if(rx.size()<1){ cerr<<"File rnd.txt did not contain valid data"<MAXRND){ cerr<<"Error: too many random multipliers ("<>mbid>>hex>>omid>>dec>>om.r[0]>>om.r[1]>>om.r[2]>>om.str>>om.dom){ om.r[2]+=zoff; i3oms.push_back(om); } inFile.close(); } } { ifstream inFile((dir+"eff-f2k").c_str(), ifstream::in); if(!inFile.fail()){ int typ; ikey om; float eff; string in; while(getline(inFile, in)){ int num=sscanf(in.c_str(), "%d %d %f %d", &om.str, &om.dom, &eff, &typ); if(num<4) typ=0; if(num>=3) if(om.isinice()) rdes.insert(make_pair(om, make_pair(eff, typ))); } inFile.close(); } } { ifstream inFile((dir+"hvs-f2k").c_str(), ifstream::in); if(!inFile.fail()){ ikey om; float hv; while(inFile>>om.str>>om.dom>>hv) if(om.isinice()) hvs.insert(make_pair(om, hv)); inFile.close(); } } { // initialize geometry vector oms; vector names; int nhqe=0; sort(i3oms.begin(), i3oms.end()); for(vector::const_iterator i=i3oms.begin(); i!=i3oms.end(); ++i) if(i->isinice()){ oms.push_back(*i); ikey om(*i); int t; float r, h; { map >::iterator j=rdes.find(om); if(j!=rdes.end()){ nhqe++; r=j->second.first; t=j->second.second; } else r=1, t=0; if(t==0){ if(r>rmax) rmax=r; } else{ if(r>rmay) rmay=r; } } if(hvs.empty()) h=1200; else{ map::iterator j=hvs.find(om); h=j==hvs.end()?0:j->second; } names.push_back(name(om, r, t, h)); } if(nhqe>0) cerr<<"Loaded "<MAXGEO){ cerr<<"Error: too many OMs ("< num; { map l; map sc; for(int n=0; nom.r[2]) l[str]=om.r[2]; if(s.hn) s.n=n; } num[str]++; s.x+=om.r[0], s.y+=om.r[1]; } if(sc.size()>NSTR){ cerr<<"Number of strings exceeds capacity of "<::iterator i=num.begin(); i!=num.end(); ++i){ unsigned char str=i->first, n=i->second; line & s = sc[str]; float d=s.h-l[str]; if(n>1 && d<=0){ cerr<<"Cannot estimate the spacing along string "<<(int)str<1?(n-1)/d:0; s.dl=0; s.dh=0; } for(int n=0; ns.r) s.r=dr; if(s.d>0){ float dz=om.r[2]-(s.h+(s.n-n)/s.d); if(s.dl>dz) s.dl=dz; if(s.dh::iterator i=num.begin(); i!=num.end(); ++i, n++){ unsigned char str=i->first; line & s = sc[str]; s.max=i->second-1; i->second=n; s.r=d.R+sqrt(s.r); if(d.rx cells[CX][CY]; { float cl[2]={0,0}, ch[2]={0,0}, crst[2]; int n=0; for(map::iterator i=num.begin(); i!=num.end(); ++i, n++){ line & s = d.sc[i->second]; for(int m=0; m<2; m++){ if(n==0 || ctr(s, m)ch[m]) ch[m]=ctr(s, m); } } d.cn[0]=CX; d.cn[1]=CY; for(int m=0; m<2; m++){ float diff=ch[m]-cl[m]; #ifdef ROMB d.cn[m]=min(d.cn[m], 1+2*(int)lroundf(diff/125)); #endif if(d.cn[m]<=1){ ch[m]=cl[m]=(cl[m]+ch[m])/2; crst[m]=1/(d.rx*(2+XXX)+diff); } else{ float s=d.R*(d.cn[m]-1); if(diff<2*s){ cerr<<"Warning: tight string packing in direction "<<(m<1?"x":"y")<::iterator i=num.begin(); i!=num.end(); ++i){ line & s = d.sc[i->second]; int n[2]; for(int m=0; m<2; m++){ n[m]=lroundf((ctr(s, m)-cl[m])*crst[m]); if(n[m]<0 && n[m]>=d.cn[m]){ cerr<<"Error in cell initialization"<first<<" too close to cell boundary"<first]++; } if(flag) d.rx=0; for(int m=0; m<2; m++){ d.cl[m]=cl[m]; d.crst[m]=crst[m]; } } { unsigned int pos=0; for(int i=0; i & c = cells[i][j]; if(c.size()>0){ d.is[i][j]=pos; for(map::const_iterator n=c.begin(); n!=c.end(); ++n){ if(pos==NSTR){ cerr<<"Number of string cells exceeds capacity of "<first]; } d.ls[pos-1]|=0x80; } else d.is[i][j]=0x80; } } cerr<<"Loaded "<> aux) d.mas=aux, n++; for(int i=0; i> aux) d.s[i]=aux, n++; if(n>0) cerr<<"Loaded "<0) d.eff*=d.mas; } #endif #ifdef TILT { d.lnum=0; d.l0=0, d.r0=0; const float cv=FPI/180, thx=225; d.lnx=cos(cv*thx), d.lny=sin(cv*thx); ifstream inFile((dir+"tilt.par").c_str(), ifstream::in); if(!inFile.fail()){ int str; float aux; vector lr; while(inFile >> str >> aux){ if(aux==0) d.l0=str; lr.push_back(aux); } inFile.close(); int size=lr.size(); if(size>LMAX){ cerr << "File tilt.par defines too many dust maps" << endl; exit(1); } for(int i=1; i pts(d.lnum), ds; vector< vector > lp(d.lnum); while(inFile >> depth){ int i=0; while(inFile >> pts[i++]) if(i>=d.lnum) break; if(i!=d.lnum) break; ds.push_back(depth); for(i=0; iLYRS){ cerr << "File tilt.dat defines too many map points" << endl; exit(1); } for(int i=1; i0) cerr<<"Loaded "< wx, wy, wz; { char * WFLA=getenv("WFLA"); float wfla=WFLA==NULL?0:atof(WFLA); if(wfla>0){ wx.push_back(0); wy.push_back(wfla-XXX); wx.push_back(1); wy.push_back(wfla+XXX); cerr<<"Using single wavelength="<>xa>>ya){ if(( xa<0 || 10 && ( 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 "< qx, qy; ifstream inFile((dir+"wv.rde").c_str(), ifstream::in); if(!inFile.fail()){ int num=0; bool flag=true; float xa, ya, yo=0; while(inFile>>ya>>xa){ if(xa<0 || (num>0 && ya<=yo)){ flag=false; break; } qx.push_back(xa); qy.push_back(ya); yo=ya; num++; } if(qx.size()<2) flag=false; inFile.close(); if(flag){ cerr<<"Loaded "<::iterator i=wy.begin(); i!=wy.end(); i++){ float w=*i, r; for(; kw) break; if(k==0) r=qx[0]; else if(k==n) r=qx[n-1]; else{ r=((w-qy[k-1])*qx[k]+(qy[k]-w)*qx[k-1])/(qy[k]-qy[k-1]); k--; } wz.push_back(r); } float eff=0, tmp=0; wx[0]=0; for(unsigned int i=1; i::iterator i=wx.begin(); i!=wx.end(); i++) *i/=eff; d.eff*=eff; } else d.eff*=max(rmax, rmay); } float wv0=400; float A, B, D, E, a, k; float Ae, Be, De, Ee, ae, ke; vector dp, be, ba, td; { bool flag=true, fail=false; ifstream inFile((dir+"icemodel.par").c_str(), ifstream::in); if((flag=!inFile.fail())){ if(flag) flag=(bool)(inFile >> a >> ae); if(flag) flag=(bool)(inFile >> k >> ke); if(flag) flag=(bool)(inFile >> A >> Ae); if(flag) flag=(bool)(inFile >> B >> Be); fail=!flag; if(flag) flag=(bool)(inFile >> D >> De); if(!flag) D=pow(wv0, k); if(flag) flag=(bool)(inFile >> E >> Ee); if(!flag) E=0; if(fail) cerr << "File icemodel.par found, but is corrupt" << endl; inFile.close(); if(fail) exit(1); } else{ cerr << "File icemodel.par was not found" << endl; exit(1); } } { ifstream inFile((dir+"icemodel.dat").c_str(), 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; i0) if(fabsf(dp[i]-dp[i-1]-dh)>dh*XXX){ cerr << "Ice table does not use uniform depth spacing" << endl; exit(1); } cerr<<"Loaded "<MAXLYS){ cerr<<"Error: too many layers ("<> bble >> bblz >> bbly)){ cerr << "File icemodel.bbl found, but is corrupt" << endl; bble=0; } else{ cerr << "Air bubble parameters: " << bble << " " << bblz << " " << bbly << endl; } } } for(int n=0; n0?dp[j]0 && abs>0) w.z[i].sca=sca, w.z[i].abs=abs; else{ 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.wvl=wva; w.ocm=ng/c; w.coschr=1/np; w.sinchr=sqrt(1-w.coschr*w.coschr); } d.fla=-1; { char * FLDR=getenv("FLDR"); d.fldr=FLDR==NULL?-1:atof(FLDR); if(d.fldr>=0){ float fold=int(d.fldr/360); float dir=d.fldr-360*fold++; cerr<<"Flasher LEDs are in a "<