#include const double pi = 3.14159; typedef double Vector[4];//Vector[0] is not used struct RotVector { double sinRot;//sine of the rotation angle double cosRot;//cosine of the rotation angle Vector axis;// the rotation axis vector }; void OrientVectorTriad(Vector &incidentVector, Vector &scatterVector, Vector &crystalAxis, double &scatAngle); double CigarEventProbability(double scatAngle,Vector crystalAxis, Vector observer, Vector scatVector, double targetTip); double NormalProbabilityDistribution(double mean, double sigma, double x); void FindRotation(Vector k1,Vector k2, RotVector &rotation); void RotateWith(RotVector omega, Vector &inVector); double Rnd();//This is an external random number generator, uniform in [0,1] double NormalProbabilityDistribution(double mean, double sigma, double x) { //This is just a gaussian with mean = mean and standard deviation = sigma return exp(-(x-mean)/sigma*(x-mean)/sigma/2.0); } void FindRotation(Vector k1,Vector k2, RotVector &rotation) { //This procedure finds the rotation vector required //to rotate unit vector k1 into unit vector k2. //First find the normalised cross product k1 X k2. That defines the rotation axis rotation.sinRot = 0; rotation.axis[1] = k1[2]*k2[3]-k1[3]*k2[2]; rotation.sinRot = rotation.sinRot + rotation.axis[1]*rotation.axis[1]; rotation.axis[2] = k1[3]*k2[1]-k1[1]*k2[3]; rotation.sinRot = rotation.sinRot + rotation.axis[2]*rotation.axis[2]; rotation.axis[3] = k1[1]*k2[2]-k1[2]*k2[1]; rotation.sinRot = rotation.sinRot + rotation.axis[3]*rotation.axis[3]; rotation.sinRot = sqrt(rotation.sinRot); //sinRot will be the sine of the rotation angle: sinRot = |k1 X k2| rotation.cosRot = 0; for (int i=1; i<=3; i++) { rotation.cosRot = rotation.cosRot + k1[i]*k2[i] ; // The scalar product k1·k2 is the cosine of the rotation angle rotation.axis[i] = rotation.axis[i]/rotation.sinRot; //Make the axis vector a unit vector } } void RotateWith(RotVector omega, Vector &inVector) { double norm = 0.0; double kw = 0.0; double e1[4]; double e2[4]; //This procedure rotates the inVector using the rotation vector omega //Find a unit vector perpendicular to omega and the vector e2[1] = omega.axis[2]*inVector[3]-omega.axis[3]*inVector[2]; norm = norm + e2[1]*e2[1]; e2[2] = omega.axis[3]*inVector[1]-omega.axis[1]*inVector[3]; norm = norm + e2[2]*e2[2]; e2[3] = omega.axis[1]*inVector[2]-omega.axis[2]*inVector[1]; norm = norm + e2[3]* e2[3]; norm = sqrt(norm); for (int i=1; i<=3; i++) e2[i] = e2[i]/norm; //Find the remaining orthogonal base vector e1[1] = e2[2]*omega.axis[3]-e2[3]*omega.axis[2]; e1[2] = e2[3]*omega.axis[1]-e2[1]*omega.axis[3]; e1[3] = e2[1]*omega.axis[2]-e2[2]*omega.axis[1]; //Here is the scalar product k·w kw = inVector[1]*omega.axis[1]+inVector[2]*omega.axis[2]+inVector[3]*omega.axis[3]; norm : sqrt(1 - kw*kw); for (int i=1; i<=3; i++) inVector[i] = kw*omega.axis[i] + e1[i]*norm*omega.cosRot + e2[i]*norm*omega.sinRot; } void OrientVectorTriad(Vector &incidentVector, Vector &scatVector, Vector &crystalAxis, double &scatAngle) //This procedure orients the event vectors in a standard way: //The incident vector points in the z direction and the //scattered vector lies in the xz plane //Note that this standard orientation differs from that defined //in the paper. Save some computing time. //For Parry simulations the crystal axis is a unit vector //perpendicular to a side face. //IMPORTANT!! //It is assumed that the raytracing preceding the entering of //this procedure has been done with tha incident ray having random //direction and with the crystal symmetry axis //perpendicular to the endfaces pointing in the z direction and //with one side face perpendicular to the x axis. { double norm; RotVector omega; Vector tempVector1; Vector tempVector2; Vector zAxis; //Set the z direction zAxis[1] = 0.0; zAxis[2] = 0.0; zAxis[3] = 1.0; //Find rotation to bring the incident ray to the positive z direction FindRotation(incidentVector, zAxis, omega); //Rotate the triad with the found rotation} RotateWith(omega, crystalAxis); RotateWith(omega, scatVector); incidentVector[1] = 0.0; incidentVector[2] = 0.0; incidentVector[3] = 1.0; //Project the scatter vector into the xy plane norm = sqrt(scatVector[1]*scatVector[1]+scatVector[2]*scatVector[2]); tempVector1[1] = scatVector[1]/norm; tempVector1[2] = scatVector[2]/norm; tempVector1[3] = 0; //This is a unit x-vector tempVector2[1] = 1.0; tempVector2[2] = 0.0; tempVector2[3] = 0.0; //Find the z axis rotation to bring the scatter vector into the xz plane FindRotation(tempVector1, tempVector2, omega); //and rotate the event vectors with the found rotation RotateWith(omega, crystalAxis); RotateWith(omega, scatVector); //The triad is now in a standard orientation and fully defined //and we are ready to enter the next procedure. //Compute the scattering angle scatAngle = 0.0; for (int i=1; i<=3; i++) scatAngle = scatAngle + incidentVector[i]*scatVector[i]; scatAngle = atan(sqrt(1-scatAngle*scatAngle)/scatAngle); if (scatAngle < 0) scatAngle = scatAngle+pi; } double CigarEventProbability(double scatAngle, Vector crystalAxis, Vector observer, Vector scatVector, double targetTip) { //This is the main procedure. Using the incident direction (0,0,1), //and the scattering angle (scatAngle)a point is generated on //the surface of the Minnaert cigar. //Then we find the rotation around the y axis necessary to make the incident //direction point at this point and rotate the crystal axis //and the scattered ray using this rotation. We then rotate //the crystal axis and the scattered ray randomly around the //axis "observer" defined by the vector from the sourceto the observer. //We then find the probability that the crystal tip now has the orientation //defined by targetTip. //This is the procedure that recycles every event a number of times //proportional to the cigar weight. double h; double theta; double alpha; double s, c, temp; double phi; double tip; double probability; Vector targetVector; RotVector omega; //Generate a latitudinal point on the cigar h = Rnd(); theta = scatAngle*h; //Find the rotation angle to bring the incident ray to this point alpha = pi/2 - atan(observer[3]/observer[1]) - theta; c = cos(alpha); s = sin(alpha); //Rotate the scattered ray and the crystal axis using this angle temp = scatVector[1]; scatVector[1] = temp*c + scatVector[3]*s; scatVector[3] = -temp*s + scatVector[3]*c; temp = crystalAxis[1]; crystalAxis[1] = temp*c + crystalAxis[3]*s; crystalAxis[3] = -temp*s + crystalAxis[3]*c; //Generate a rotation angle for an axial rotation on the cigar phi = 2*pi* Rnd(); //Set up a rotation vector. The rotation axis points to the observer omega.cosRot = cos(phi); omega.sinRot = sin(phi); omega.axis[1] = observer[1]; omega.axis[2] = observer[2]; omega.axis[3] = observer[3]; //Rotate the scattered ray and the crystal axis using the rotation vector RotateWith(omega, scatVector); RotateWith(omega, crystalAxis); //Find the tip angle} tip = fabs(crystalAxis[3]); tip = atan(sqrt(1-tip*tip)/tip); //Check the tip probability distribution and generate the event probability. //Weight the event probability with the factor scatAngle/pi. //Target tip angles below 40 degrees are here considered to mean plate crystals, //angles above 50 degrees mean column crystals and in the interval //40 to 50 degrees mean randomly oriented crystals. //Parry simulations may require a somewhat more elaborate code if (targetTip < 40*pi/180) probability = NormalProbabilityDistribution(0, targetTip, tip); //Plates else if (targetTip > 50*pi/180) probability = NormalProbabilityDistribution(pi/2, pi/2-targetTip, tip); //Columns else probability = 1.0; //Random crystals return probability; }