#include <cmath>
#include <cstdlib>
#include <iostream>
#include <map>
#include <vector>
#include <fstream>

using namespace std;

#include "ice.cxx"

double xrnd(){
  double rnd=0;
  while(rnd==0){
    rnd=rand();
    rnd/=RAND_MAX;
  }
  return rnd;
}

struct KEY{
  int str, om;

  KEY(){
    str=-1; om=-1;
  }

  bool operator< (const KEY& rhs) const{
    return str < rhs.str || (str == rhs.str && om < rhs.om);
  }
};

struct STR{
  bool zf;
  int l, h;
  double zn, zx, sp;
  double rl[3], rh[3];

  STR(){
    zf=true;
    sp=-1;
  }
};

map<int,STR> strs;

struct DOM{
  double r[3], R;

  DOM(){
    R=0.178;
    R*=5;
  }
};
map<KEY,DOM> doms;

void load(){
  ifstream inFile;
  inFile.open("geo-f2k", ifstream::in);

  if(!inFile.fail()){
    KEY key;
    DOM dom;
    string id;
    unsigned long long domid;
    while(inFile>>id>>hex>>domid>>dec>>dom.r[0]>>dom.r[1]>>dom.r[2]>>key.str>>key.om) if(key.str>0 && key.om>=1 && key.om<=60){
      doms[key]=dom;
      STR & s=strs[key.str];
      if(s.sp<0){
	s.l=key.om; s.h=key.om;
	for(int i=0; i<3; i++){
	  s.rl[i]=dom.r[i];
	  s.rh[i]=dom.r[i];
	}
	s.sp=0;
      }
      else{
	if(key.om<s.l) s.l=key.om;
	if(key.om>s.h) s.h=key.om;
	for(int i=0; i<3; i++){
	  if(dom.r[i]<s.rl[i]) s.rl[i]=dom.r[i];
	  if(dom.r[i]>s.rh[i]) s.rh[i]=dom.r[i];
	}
	s.sp++;
      }
    }
    inFile.close();
  }

  for(map<KEY,DOM>::iterator s=doms.begin(); s!=doms.end(); s++){
    STR & str=strs[s->first.str];

    double dz=s->second.r[2]-(str.rl[2]+(s->first.om-str.h)*(str.rh[2]-str.rl[2])/(str.l-str.h));
    if(str.zf){
      str.zn=dz;
      str.zx=dz;
      str.zf=false;
    }
    else{
      if(dz<str.zn) str.zn=dz;
      if(dz>str.zx) str.zx=dz;
    }
  }

  cerr<<doms.size()<<" DOMs loaded on "<<strs.size()<<" strings"<<endl;
}

struct SRC{
  double r[3], n[3];

  SRC(){
    load();
    for(int i=0; i<3; i++){ r[i]=i<2?0:-1950; n[i]=i==0?1:0; }
  }

  void init(int str, int om){
    for(map<KEY,DOM>::iterator s=doms.begin(); s!=doms.end(); s++) if(s->first.str==str && s->first.om==om){
      for(int i=0; i<3; i++) r[i]=s->second.r[i]; doms.erase(s); break;
    }    
  }

  void circ(){
    static const double rms=9.2*sqrt(3)*M_PI/180;
    double xi=xrnd();
    xi=rms*(2*xi-1);
    double nz=sin(xi), np=cos(xi);

    xi=xrnd();
    xi*=2*M_PI;
    n[0]=np*cos(xi);
    n[1]=np*sin(xi);
    n[2]=nz;
  }
} x;

struct quat{
  double val[4];
};

quat operator*(const quat &l, const quat &r){
  quat z;
  const double & l0=l.val[0], l1=l.val[1], l2=l.val[2], l3=l.val[3];
  const double & r0=r.val[0], r1=r.val[1], r2=r.val[2], r3=r.val[3];
  z.val[0]=l0*r0-l1*r1-l2*r2-l3*r3;
  z.val[1]=l0*r1+l1*r0+l2*r3-l3*r2;
  z.val[2]=l0*r2-l1*r3+l2*r0+l3*r1;
  z.val[3]=l0*r3+l1*r2-l2*r1+l3*r0;
  return z;
}

struct photon{
  double r[3], n[3], d, t;
  double c, g, la, ls;
  double xx;

  photon(){
    d=0; t=0;
    for(int i=0; i<3; i++){
      r[i]=x.r[i]; n[i]=x.n[i];
    }

    xx=1.e-5;
    c=m.c; g=m.g;
    la=m.la; ls=m.ls;
  }

  void rotate(){
    double u[3];
    {
      double p1[3], p2[3];
      {
	int i0=0;
	double n0=n[0];
	for(int i=1; i<3; i++) if(abs(n[i])>abs(n0)){ i0=i; n0=n[i]; }
	int i1=(i0+1)%3, i2=(i0+2)%3;
	double n1=n[i1], n2=n[i2], nq=n0*n0;
	double r1=sqrt(nq+n1*n1), r2=sqrt(nq+n2*n2);

	p1[i0]=-n1/r1; p1[i1]=n0/r1; p1[i2]=0;
	p2[i0]=-n2/r2; p2[i1]=0; p2[i2]=n0/r2;
      }
      {
	double q1[3], q2[3], r1=0, r2=0;
	for(int i=0; i<3; i++){
	  double a1=p1[i]-p2[i]; q1[i]=a1; r1+=a1*a1;
	  double a2=p1[i]+p2[i]; q2[i]=a2; r2+=a2*a2;
	}

	r1=sqrt(r1); r2=sqrt(r2);

	for(int i=0; i<3; i++){
	  p1[i]=q1[i]/r1;
	  p2[i]=q2[i]/r2;
	}
      }

      {
	double zi=xrnd();
	zi*=2*M_PI;
	double px=cos(zi), py=sin(zi);
	for(int i=0; i<3; i++) u[i]=px*p1[i]+py*p2[i];
      }
    }

    {
      double xi=xrnd();
      xi=2*xi-1;
      if(g!=0){
	double g2=g*g;
	double ga=(1-g2)/(1+g*xi);
	xi=(1+g2-ga*ga)/(2*g);
      }
      double c=sqrt((1+xi)/2), s=sqrt((1-xi)/2);
      quat q1, q2, v;
      {
	q1.val[0]=c;
	q2.val[0]=c;
	v.val[0]=0;
      }
      for(int i=0; i<3; i++){
	double a=s*u[i];
	q1.val[i+1]=a;
	q2.val[i+1]=-a;
	v.val[i+1]=n[i];
      }
      quat f=q1*v*q2;
      double r=0;
      for(int i=0; i<3; i++){
	double a=f.val[i+1];
	n[i]=a; r+=a*a;
      }
      if(abs(1-r)>xx){
	r=sqrt(r);
	for(int i=0; i<3; i++) n[i]/=r;
      }
    }
  }

  void advance(double a){
    d+=a; t+=a/c;
    for(int i=0; i<3; i++) r[i]+=n[i]*a;
  }

  void propagate(){
    double TOT=-log(xrnd());
    while(true){
      double tot=m.aX(TOT, r[2], n[2]);
      double SCA=-log(xrnd());
      double sca=m.sX(SCA, r[2], n[2]);
      if(sca>tot) sca=tot;

      KEY om=sphere(sca);
      TOT-=m.xA(sca, r[2], n[2]);
      advance(sca);

      if(om.str>=0){
	cout<<om.str<<" "<<om.om<<"  "<<t<<" "<<n[2]<<endl;
	break;
      }

      if(abs(TOT)<xx) TOT=0;
      if(TOT==0) break;

      rotate();
    }
  }

  KEY sphere(double &sca){
    KEY imin;
    double hmin;
    for(map<int,STR>::const_iterator str=strs.begin(); str!=strs.end(); str++){
      bool pass=false;
      const STR & s=str->second;
      for(int i=0; i<3; i++){
	if(n[i]>0){
	  if(s.rh[i]<r[i] || s.rl[i]>r[i]+sca*n[i]){
	    pass=true;
	    break;
	  }
	}
	else{
	  if(s.rl[i]>r[i] || s.rh[i]<r[i]+sca*n[i]){
	    pass=true;
	    break;
	  }
	}
      }
      if(pass) continue;

      double min, max;
      double st=(s.h-s.l)/(s.rh[2]-s.rl[2]);
      if(n[2]>0){
	min=s.h-(r[2]+sca*n[2]-s.zn-s.rl[2])*st;
	max=s.h-(r[2]-s.zx-s.rl[2])*st;
      }
      else{
	min=s.h-(r[2]-s.zn-s.rl[2])*st;
	max=s.h-(r[2]+sca*n[2]-s.zx-s.rl[2])*st;
      }
      int n1=(int)floor(min);
      int n2=(int)ceil(max);
      if(n1>s.h | n2<s.l) continue;

      for(int om=n1; om<=n2; om++){
	KEY key;
	key.str=str->first;
	key.om=om;
	DOM dom=doms[key];

	double R=dom.R;
	double b=0, c=0;
	for(int i=0; i<3; i++){
	  double dr=dom.r[i]-r[i];
	  b+=n[i]*dr; c+=dr*dr;
	}
	c-=R*R;
	double D=b*b-c, h=-1;
	if(D>0){
	  D=sqrt(D);
	  h=b-D; if(h<0) h=b+D;
	}

	if(h>0 && h<sca) if(imin.str<0 || h<hmin){ imin=key; hmin=h; }
      }
    }

    if(imin.str>=0) sca=hmin;
    return imin;
  }
};

main(int arg_c, char *arg_a[]){
  if(arg_c>1){
    int str=0, dom=0;
    unsigned int seed=0;
    unsigned long long num=1000000ULL;

    if(arg_c>1) str=atoi(arg_a[1]);
    if(arg_c>2) dom=atoi(arg_a[2]);
    if(arg_c>3) num=(unsigned long long) atof(arg_a[3]);
    if(arg_c>4) seed=atoi(arg_a[4]);

    srand(1+seed);
    x.init(str, dom);
    for(unsigned long long i=0; i<num; i++){
      x.circ();
      photon p;
      p.propagate();
    }
  }
  else cerr<<"Use: "<<arg_a[0]<<" [str] [om] [num] [seed]"<<endl;
}