Adjusted pulse profile computation to reflect new model

* Using variables as defined in Ben's formulas
* Using axes relations and rotation directions as defined in Ben's formulas

src/pulsaranimationwidget.cpp

@@ -20,6 +20,9 @@
 
 #include "pulsaranimationwidget.h"
 
+// workaround for lack of C99 support in C++
+#define isnan(x) ((x) != (x))
+
 #define GL_PI 3.1415f
 const double PulsarAnimationWidget::deg2rad = PI/180.0;
 
@@ -332,7 +335,7 @@ void PulsarAnimationWidget::paintGL()
             glRotatef(-m_pulsarSpinAxisInclination, 1.0, 0.0, 0.0);
 
             glRotatef(m_pulsarRotationAngle - 90.0, 0.0, 0.0, 1.0);
-            glRotatef(-m_pulsarMagneticAxisInclination, 0.0, 1.0, 0.0);
+            glRotatef(m_pulsarMagneticAxisInclination, 0.0, 1.0, 0.0);
 
             // draw magnetic axis (for both cones)
             if(m_showRotationAxes) {
@@ -412,7 +415,7 @@ void PulsarAnimationWidget::paintGL()
             glRotatef(-m_pulsarSpinAxisInclination, 1.0, 0.0, 0.0);
 
             glRotatef(-m_pulsarRotationAngle - 90.0, 0.0, 0.0, 1.0);
-            glRotatef(-m_pulsarMagneticAxisInclination, 0.0, 1.0, 0.0);
+            glRotatef(m_pulsarMagneticAxisInclination, 0.0, 1.0, 0.0);
 
             glTranslatef(0.0, 0.0, -4.0);
 
@@ -603,7 +606,7 @@ void PulsarAnimationWidget::updateCameraPosition(const double angleH, const doub
 
     updateGL();
 
-    //	qDebug("Camera (x,y,z,+,h,v): %f, %f, %f, %f, %f, %f", m_cameraPosX, m_cameraPosY, m_cameraPosZ, zoom, angleH, angleV);
+    qDebug("Camera (x,y,z,+,h,v): %f, %f, %f, %f, %f, %f", m_cameraPosX, m_cameraPosY, m_cameraPosZ, zoom, angleH, angleV);
 }
 
 void PulsarAnimationWidget::setFramePerSecond(const unsigned int fps)
@@ -635,6 +638,8 @@ void PulsarAnimationWidget::setPulsarMagneticAxisInclination(const int degrees)
 
 void PulsarAnimationWidget::setPulsarBeamAngle(const int degrees)
 {
+    m_pulsarBeamAngle = degrees;
+
     m_pulsarBeamOuterRadius = tan(deg2rad * degrees * 0.5f) * m_pulsarBeamLength;
     m_pulsarBeamInnerRadius = m_pulsarBeamOuterRadius - m_pulsarBeamRimSize;
     if(m_pulsarBeamInnerRadius < 0.0) m_pulsarBeamInnerRadius = 0.0;
@@ -654,33 +659,30 @@ void PulsarAnimationWidget::resetParameters()
 void PulsarAnimationWidget::updatePulseProfile()
 {
     // prepare parameters (e.g. convert to radians where necessary)
-    const double    i               = deg2rad * m_pulsarSpinAxisInclination;
-    const double    y               = deg2rad * m_pulsarMagneticAxisInclination;
-    double          phiOrb          = 0.0;
     const double    deltaPhiRot     = deg2rad * 1.0;
-    const double    deltaPhiOrb     = deg2rad * deltaPhiRot;
-    const double    rp              = 0.0;
-    const double    xk              = -m_cameraPosZ;
-    const double    yk              = -m_cameraPosX;
-    const double    zk              = m_cameraPosY;
-    const double    cam             = pow(xk, 2.0) + pow(yk, 2.0) + pow(zk, 2.0);
-    const double    alpha           = deg2rad * (90.0 - m_mouseAngleH);
-    const double    delta           = deg2rad * m_mouseAngleV;
-    const double    gaussProfile    = 0.012337;
+    const double    gamma           = deg2rad * m_pulsarSpinAxisInclination;
+    const double    alpha           = deg2rad * m_pulsarMagneticAxisInclination;
+    const double    beta            = alpha - gamma;
+    const double    rho             = deg2rad * m_pulsarBeamAngle * 0.5;
+    const double    gaussProfile    = 0.04;
 
-    for(int x = 0; x < 360; ++x) {
-        // determine angle between pulsar's magnetic axis and line of sight
-        phiOrb += deltaPhiOrb;
-        const double phiRot = x * deltaPhiRot;
+    const double    t1              = pow(sin(rho*0.5), 2) - pow(sin(beta*0.5), 2);
+    const double    t2              = sin(alpha) * sin(beta+alpha);
+    const double    W               = 2 * asin(sqrt( t1 / t2 ));
 
-        double a = -sin(y) * sin(phiRot) * (xk + rp * cos(phiOrb)) \
-                   + (cos(i) * sin(y) * cos(phiRot) + sin(i) * cos(y)) * (yk + rp * sin(phiOrb)) \
-                   - (sin(i) * sin(y) * cos(phiRot) - cos(i) * cos(y)) * zk;
+    qDebug("alpha: %f, beta: %f, rho: %f, t1: %f, t2: %f, W: %f", alpha*180/PI, beta*180/PI, rho*180/PI, t1*180/PI, t2*180/PI, W*180/PI);
 
-        double b = sqrt(pow(rp,2.0) + cam - (2.0 * sqrt(cam) * rp * cos(delta) * sin(alpha + phiOrb)));
+    for(int x = 0; x < 360; ++x) {
+        // update pulsar rotation angle
+        const double phiRot = x * deltaPhiRot;
 
         // determine and store pulse amplitude
-        m_pulseProfile[x] = exp(-2.0 * (1.0 - fabs(a/b)) / gaussProfile);
+        if(!isnan(W)){
+            m_pulseProfile[x] = exp(-2.0 * pow(phiRot -W, 2) / gaussProfile) + exp(-2.0 * pow(phiRot + W - 2.0 * PI, 2) / gaussProfile);
+        }
+        else {
+            m_pulseProfile[x] = 0.0;
+        }
     }
 
     // propagate new profile

src/pulsaranimationwidget.h

@@ -96,6 +96,7 @@ private:
     double m_pulsarSpinAxisInclination;
     double m_pulsarMagneticAxisInclination;
     double m_pulsarBeamLength;
+    double m_pulsarBeamAngle;
     double m_pulsarBeamInnerRadius;
     double m_pulsarBeamOuterRadius;
     double m_pulsarBeamRimSize;