Index: trunk/Mars/hawc/transition.C
===================================================================
--- trunk/Mars/hawc/transition.C	(revision 19631)
+++ trunk/Mars/hawc/transition.C	(revision 19631)
@@ -0,0 +1,108 @@
+Bool_t HandleInput()
+{
+    TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
+    while (1)
+    {
+        //
+        // While reading the input process gui events asynchronously
+        //
+        timer.TurnOn();
+        TString input = Getline("Type 'q' to exit, <return> to go on: ");
+        timer.TurnOff();
+
+        if (input=="q\n")
+            return kFALSE;
+
+        if (input=="\n")
+            return kTRUE;
+    };
+
+    return kFALSE;
+}
+
+void PropagateZ(TVector3 &p, TVector3 &w, double z, TArrow &arr)
+{
+    arr.SetX1(p.X());
+    arr.SetY1(p.Z());
+
+    p += (z-p.Z())/w.Z()*w;
+
+    arr.SetX2(p.X());
+    arr.SetY2(p.Z());
+}
+
+void PropagateDZ(TVector3 &p, TVector3 &w, double dz, TArrow &arr)
+{
+    arr.SetX1(p.X());
+    arr.SetY1(p.Z());
+
+    p += dz/w.Z()*w;
+
+    arr.SetX2(p.X());
+    arr.SetY2(p.Z());
+}
+
+void transition()
+{
+    TH1C h("", "", 100, -2, 2);
+    h.SetMinimum(-1.5);
+    h.SetMaximum(1.5);
+    h.SetStats(kFALSE);
+    h.DrawCopy();
+
+    TLine line;
+    line.DrawLine(-2, 0, 2, 0);
+
+    TVector3 norm(0, 0, -1);
+    norm *= 1./norm.Mag();
+
+    TArrow arr;
+    arr.SetArrowSize(0.01);
+    arr.DrawArrow(0, 0, 0+norm.X(), 0+norm.Z());
+
+    TArrow arr_in(1, 1, 0, 0);
+    arr_in.SetArrowSize(0.01);
+    arr_in.SetLineColor(kBlue);
+    arr_in.Draw();
+
+    TArrow arr_out(0, 0, 1, 1);
+    arr_out.SetArrowSize(0.01);
+    arr_out.SetLineColor(kRed);
+    arr_out.Draw();
+
+    while (1)
+    {
+        double x, y;
+        gRandom->Circle(x, y, 1);
+
+        TVector3 pos(x, 0, y);
+        TVector3 dir = -pos;
+
+        arr_in.SetX1(-dir.X()*1.5);
+        arr_in.SetY1(-dir.Z()*1.5);
+
+        double n1 = dir.Z()>0 ? 1.5 : 1;   // where the ray comes from
+        double n2 = dir.Z()>0 ? 1   : 1.5; // where the ray goes to
+
+        cout << sin(acos(dir*norm)) << endl;
+
+        int rc = MOptics::ApplyTransition(dir, norm, n1, n2);
+
+        cout << rc << endl;
+
+        if (rc>=3) // Refraction
+        {
+            // Snell's law: n2*sin(theta2) = n1*sin(theta1) -- Is thge result correct?
+            cout << n2/n1*sin(acos(dir*norm)) << endl;
+        }
+
+        arr_out.SetX2(dir.X()*1.5);
+        arr_out.SetY2(dir.Z()*1.5);
+
+        gPad->Modified();
+        gPad->Update();
+
+        if (!HandleInput())
+            break;
+    }
+}