Index: /trunk/Mars/msimreflector/MSimRays.cc
===================================================================
--- /trunk/Mars/msimreflector/MSimRays.cc	(revision 19787)
+++ /trunk/Mars/msimreflector/MSimRays.cc	(revision 19788)
@@ -165,21 +165,46 @@
 
         Double_t x, y;
-        const Double_t r = gRandom->Uniform();
-        gRandom->Circle(x, y, maxr*TMath::Sqrt(r));
-/*
-        Double_t ra = gRandom->Uniform(maxr);
-        Double_t ph = gRandom->Uniform(TMath::TwoPi());
-
-
-        // Get radom incident point on the mirror plane.
-        //const Double_t x = gRandom->Uniform(-maxr, maxr);
-        //const Double_t y = gRandom->Uniform(-maxr, maxr);
-
-        Double_t x = ra*sin(ph);
-        Double_t y = ra*cos(ph);
-
-//        if (x*x + y*y > maxr*maxr)
-//            continue;
-  */
+        if (fHeight<0)
+        {
+            // Parallel light
+            // --------------
+            const Double_t r = gRandom->Uniform();
+            gRandom->Circle(x, y, maxr*TMath::Sqrt(r));
+        }
+        else
+        {
+            // Point source
+            // ------------
+            // Adapted from: http://mathworld.wolfram.com/SpherePointPicking.html
+            // Note that theta and phi is exchanged!
+ 
+            // The maximum zenith angle is theta=atan(maxr/h)
+            // cos(theta) = cos(atan(maxr/h)) = 1/sqrt(1+maxr^2/h^2)
+            const double min_cost  = 1./TMath::Sqrt(1.+maxr*maxr/h/h);
+            const double cos_theta = gRandom->Uniform(min_cost, 1);
+
+            gRandom->Circle(x, y, h*TMath::Sqrt(1./cos_theta/cos_theta - 1));
+
+            // const double cos_theta = gRandom->Uniform(ct, 1);
+            // const double sin_theta = TMath::Sqrt(1.-cos_theta*cos_theta);
+            // gRandom->Circle(x, y, h*sin_theta/cos_theta);
+
+            // Homogeneous on a sphere
+            // const double phi = TMath::TwoPi() * gRandom->Uniform();
+            // x = sin_theta * cos(phi);
+            // y = sin_theta * sin(phi);
+            // z = cos_theta;
+
+            // Project the photons to a plane at z=1
+            // x /= cos_theta;
+            // y /= cos_theta;
+            // z /= cos_theta; // z = 1
+
+            // The radius of the sphere is h
+            // x *= h;
+            // y *= h;
+            // z *= h; // z = h
+        }
+
         // The is the incident direction of the photon
         // h==0 means infinitiy