New rotation

Makefile

@@ -49,6 +49,7 @@ endif
 ##
 
 TARGET  = lime.x 
+#CC		= gcc -fopenmp -g
 CC		= gcc -fopenmp
 SRCS	= src/aux.c src/messages.c src/grid.c src/LTEsolution.c	\
 	  src/main.c src/molinit.c src/photon.c src/popsin.c	\

doc/usermanual.rst

@@ -515,15 +515,6 @@ number.
 The image resolution or size of each pixel. This number is given in arc
 seconds. The image field of view is therefore pxls x imgres.
 
-.. code:: c
-
-    (double) img[i]->theta (required)
-
-Theta is the viewing angle (the angle between the model z axis and the
-ray-tracers line of sight). This number is given in radians, not
-degrees, so that a face-on view (of models where this term is
-applicable) is 0 and edge-on view is π/2. Note that you can use the predefined PI macro: e.g. to express π/2, write PI/2.0 in your model file.
-
 .. code:: c
 
     (double) img[i]->distance (required)
@@ -547,16 +538,6 @@ an optical depth image cube (dimensionless).
 
 This variable is the name of the output FITS file. It has no default value.
 
-.. code:: c
-
-    (double) img[i]->phi (optional)
-
-Phi is an optional geometric parameter. Like theta, it should be given
-in radians between 0 and 2π. Phi is the rotation angle of the model
-(x,y)-plane around the z-axis. If the model is view face-on (so that the
-line of sight coincides with the z-axis), phi corresponds to the
-position angle on the sky. The default value is 0.
-
 .. code:: c
 
     (double) img[i]->source_vel (optional)
@@ -633,18 +614,68 @@ bandwidth should be set. Any other combination will produce an error.
 For line images, at least one moldatfile should be provided and
 optionally a dust opacity table as well.
 
+Image rotation parameters
+~~~~~~~~~~~~~~~~~~~~~~~~~
+
+There are now two ways to specify the desired orientation of the model at the raytracing step: we have retained the old theta/phi angles, but have now added a new triplet: azimuth/inclination/PA. None of these five parameters is now mandatory. If none are provided, theta=phi=0 will be assumed. If you provide all three azimuth/inclination/PA values, these will be used instead of theta/phi, regardless if you also set either or both of theta/phi.
+
+Note that all of these angles should be given in radians. You can however use the predefined PI macro for this: e.g. to express π/2, write PI/2.0 in your model file.
+
+The rotation parameters in detail:
+
+.. code:: c
+
+    (double) img[i]->theta (optional)
+
+Theta is the vertical viewing angle (the vertical angle between the model z axis and the
+ray-tracer's line of sight). A face-on view (of models where this term is
+applicable) is 0 and edge-on view is π/2. The default value is 0.
+
+.. code:: c
+
+    (double) img[i]->phi (optional)
+
+Phi is the horizontal viewing angle (the horizontal angle between the model z axis and the
+ray-tracer's line of sight). A face-on view (of models where this term is
+applicable) is 0 and edge-on view is π/2. The default value is 0.
+
+If theta/phi are applied, for zero values of both the model X axis points to the left, Y points upward and Z points in the direction of gaze of the observer (i.e. away from the observer).
+
+.. code:: c
+
+    (double) img[i]->azimuth (optional)
+
+Azimuth rotates the model from Y towards X.
+
+.. code:: c
+
+    (double) img[i]->incl (optional)
+
+Inclination rotates the model from Z towards X.
+
+.. code:: c
+
+    (double) img[i]->posang (optional)
+
+Position angle rotates the model from Y towards X.
+
+If azimuth/incl/posang are applied (i.e. if all three values are supplied in your model.c), for zero values of all the model X axis points downward, Y points toward the right and Z towards the observer.
+
+
 Model functions
 ---------------
 
 The second part of the model.c file contains the actual model
-description. This is provided as seven subroutines: density, molecular
+description. This is provided as eight subroutines: density, molecular
 abundance, temperature, systematic velocities, random velocities,
-magnetic field, and gas-to-dust ratio. The user only needs to provide
+magnetic field, gas-to-dust ratio, and grid-point number density. The user only needs to provide
 the functions that are relevant to a particular model, e.g., for
 continuum images only, the user need not include the abundance function
 or any of the velocity functions. The magnetic field function needs only
 be included for continuum polarization images.
 
+Note that you should avoid singularities in these functions - i.e., places where LIME might attempt to divide by zero, or in some other way generate an overflow.
+
 
 .. figure:: images/fig_coords_big.png
    :alt: coordinates

example/model.c

@@ -90,11 +90,13 @@ input(inputPars *par, image *img){
   img[i].trans                  = 3;              // zero-indexed J quantum number
   img[i].pxls                   = 100;            // Pixels per dimension
   img[i].imgres                 = 0.1;            // Resolution in arc seconds
-  img[i].theta                  = 0.0;            // 0: face-on, pi/2: edge-on
   img[i].distance               = 140*PC;         // source distance in m
   img[i].source_vel             = 0;              // source velocity in m/s
   img[i].unit                   = 0;              // 0:Kelvin 1:Jansky/pixel 2:SI 3:Lsun/pixel 4:tau
   img[i].filename               = "image0.fits";  // Output filename
+  img[i].azimuth                = 0.0;
+  img[i].incl                   = 0.0;
+  img[i].posang                 = 0.0;
 }
 
 /******************************************************************************/

src/aux.c

@@ -17,10 +17,14 @@
 
 void
 parseInput(inputPars inpar, configInfo *par, image **img, molData **m){
-  int i,id;
+  int i,j,id;
   double BB[3],normBSquared,dens[MAX_N_COLL_PART];
-  double cosPhi,sinPhi,cosTheta,sinTheta,dummyVel[DIM];
+  double dummyVel[DIM];
   FILE *fp;
+  _Bool doThetaPhi;
+  double cos_pa,sin_pa,cosPhi,sinPhi,cos_incl,sin_incl,cosTheta,sinTheta,cos_az,sin_az;
+  double tempRotMat[3][3],auxRotMat[3][3];
+  int row,col;
 
   /* Copy over user-set parameters to the configInfo versions. (This seems like duplicated effort but it is a good principle to separate the two structs, for several reasons, as follows. (i) We will usually want more config parameters than user-settable ones. The separation leaves it clearer which things the user needs to (or can) set. (ii) The separation allows checking and screening out of impossible combinations of parameters. (iii) We can adopt new names (for clarity) for config parameters without bothering the user with a changed interface.) */
   par->radius       = inpar.radius;
@@ -222,38 +226,144 @@ The presence of one of these combinations at least is checked here, although the
       (*img)[i].pixel[id].tau = malloc(sizeof(double)*(*img)[i].nchan);
     }
 
-    /* Rotation matrix
+    /*
+The image rotation matrix is used within traceray() to transform the coordinates of a vector (actually two vectors - the ray direction and its starting point) as initially specified in the observer-anchored frame into the coordinate frame of the model. In linear algebra terms, the model-frame vector v_mod is related to the vector v_obs as expressed in observer- (or image-) frame coordinates via the image rotation matrix R by
+
+    v_mod = R * v_obs,				1
+
+the multiplication being the usual matrix-vector style. Note that the ith row of R is the ith axis of the model frame with coordinate values expressed in terms of the observer frame.
+
+The matrix R can be broken into a sequence of several (3 at least are needed for full degrees of freedom) simpler rotations. Since these constituent rotations are usually easier to conceive in terms of rotations of the model in the observer framework, it is convenient to invert equation (1) to give
+
+    v_obs = R^T * v_mod,			2
+
+where ^T here denotes transpose. Supposing now we rotate the model in a sequence R_3^T followed by R_2^T followed by R_1^T, equation (2) can be expanded to give
+
+    v_obs = R_1^T * R_2^T * R_3^T * v_mod.	3
+
+Inverting everything to return to the format of equation (1), which is what we need, we find
 
-            |1          0           0   |
-     R_x(a)=|0        cos(a)      sin(a)|
-            |0       -sin(a)      cos(a)|
+    v_mod = R_3 * R_2 * R_1 * v_obs.		4
 
-            |cos(b)     0       -sin(b)|
-     R_y(b)=|  0        1          0   |
-            |sin(b)     0        cos(b)|
+LIME provides two different schemes of {R_1, R_2, R_3}: {PA, phi, theta} and {PA, inclination, azimuth}. As an example, consider phi, which is a rotation of the model from the observer Z axis towards the X. The matching obs->mod rotation matrix is therefore
 
-            |      cos(b)       0          -sin(b)|
-     Rot =  |sin(a)sin(b)     cos(a)  sin(a)cos(b)|
-            |cos(a)sin(b)    -sin(a)  cos(a)cos(b)|
+            ( cos(ph)  0  -sin(ph) )
+            (                      )
+    R_phi = (    0     0     1     ).
+            (                      )
+            ( sin(ph)  0   cos(ph) )
 
     */
 
-    cosPhi   = cos((*img)[i].phi);
-    sinPhi   = sin((*img)[i].phi);
-    cosTheta = cos((*img)[i].theta);
-    sinTheta = sin((*img)[i].theta);
-    (*img)[i].rotMat[0][0] =           cosPhi;
-    (*img)[i].rotMat[0][1] =  0.0;
-    (*img)[i].rotMat[0][2] =          -sinPhi;
-    (*img)[i].rotMat[1][0] =  sinTheta*sinPhi;
-    (*img)[i].rotMat[1][1] =  cosTheta;
-    (*img)[i].rotMat[1][2] =  sinTheta*cosPhi;
-    (*img)[i].rotMat[2][0] =  cosTheta*sinPhi;
-    (*img)[i].rotMat[2][1] = -sinTheta;
-    (*img)[i].rotMat[2][2] =  cosTheta*cosPhi;
+    doThetaPhi = (((*img)[i].incl<-900.)||((*img)[i].azimuth<-900.)||((*img)[i].posang<-900.))?1:0;
+
+    if(doThetaPhi){
+      /* For the present PA is not implemented for the theta/phi scheme. Thus we just load the identity matrix at present.
+      */
+      (*img)[i].rotMat[0][0] = 1.0;
+      (*img)[i].rotMat[0][1] = 0.0;
+      (*img)[i].rotMat[0][2] = 0.0;
+      (*img)[i].rotMat[1][0] = 0.0;
+      (*img)[i].rotMat[1][1] = 1.0;
+      (*img)[i].rotMat[1][2] = 0.0;
+      (*img)[i].rotMat[2][0] = 0.0;
+      (*img)[i].rotMat[2][1] = 0.0;
+      (*img)[i].rotMat[2][2] = 1.0;
+    }else{
+      /* Load PA rotation matrix R_PA:
+      */
+      cos_pa   = cos((*img)[i].posang);
+      sin_pa   = sin((*img)[i].posang);
+      (*img)[i].rotMat[0][0] =  cos_pa;
+      (*img)[i].rotMat[0][1] = -sin_pa;
+      (*img)[i].rotMat[0][2] =  0.0;
+      (*img)[i].rotMat[1][0] =  sin_pa;
+      (*img)[i].rotMat[1][1] =  cos_pa;
+      (*img)[i].rotMat[1][2] =  0.0;
+      (*img)[i].rotMat[2][0] =  0.0;
+      (*img)[i].rotMat[2][1] =  0.0;
+      (*img)[i].rotMat[2][2] =  1.0;
+    }
+
+    if(doThetaPhi){
+      /* Load phi rotation matrix R_phi:
+      */
+      cosPhi   = cos((*img)[i].phi);
+      sinPhi   = sin((*img)[i].phi);
+      auxRotMat[0][0] =  cosPhi;
+      auxRotMat[0][1] =  0.0;
+      auxRotMat[0][2] = -sinPhi;
+      auxRotMat[1][0] =  0.0;
+      auxRotMat[1][1] =  1.0;
+      auxRotMat[1][2] =  0.0;
+      auxRotMat[2][0] =  sinPhi;
+      auxRotMat[2][1] =  0.0;
+      auxRotMat[2][2] =  cosPhi;
+    }else{
+      /* Load inclination rotation matrix R_inc:
+      */
+      cos_incl = cos((*img)[i].incl + PI);
+      sin_incl = sin((*img)[i].incl + PI);
+      auxRotMat[0][0] =  cos_incl;
+      auxRotMat[0][1] =  0.0;
+      auxRotMat[0][2] = -sin_incl;
+      auxRotMat[1][0] =  0.0;
+      auxRotMat[1][1] =  1.0;
+      auxRotMat[1][2] =  0.0;
+      auxRotMat[2][0] =  sin_incl;
+      auxRotMat[2][1] =  0.0;
+      auxRotMat[2][2] =  cos_incl;
+    }
+
+    for(row=0;row<3;row++){
+      for(col=0;col<3;col++){
+        tempRotMat[row][col] = 0.0;
+        for(j=0;j<3;j++)
+          tempRotMat[row][col] += auxRotMat[row][j]*(*img)[i].rotMat[j][col];
+      }
+    }
+
+    if(doThetaPhi){
+      /* Load theta rotation matrix R_theta:
+      */
+      cosTheta = cos((*img)[i].theta);
+      sinTheta = sin((*img)[i].theta);
+      auxRotMat[0][0] =  1.0;
+      auxRotMat[0][1] =  0.0;
+      auxRotMat[0][2] =  0.0;
+      auxRotMat[1][0] =  0.0;
+      auxRotMat[1][1] =  cosTheta;
+      auxRotMat[1][2] =  sinTheta;
+      auxRotMat[2][0] =  0.0;
+      auxRotMat[2][1] = -sinTheta;
+      auxRotMat[2][2] =  cosTheta;
+    }else{
+      /* Load azimuth rotation matrix R_az:
+      */
+      cos_az   = cos((*img)[i].azimuth + PI/2.0);
+      sin_az   = sin((*img)[i].azimuth + PI/2.0);
+      auxRotMat[0][0] =  cos_az;
+      auxRotMat[0][1] = -sin_az;
+      auxRotMat[0][2] =  0.0;
+      auxRotMat[1][0] =  sin_az;
+      auxRotMat[1][1] =  cos_az;
+      auxRotMat[1][2] =  0.0;
+      auxRotMat[2][0] =  0.0;
+      auxRotMat[2][1] =  0.0;
+      auxRotMat[2][2] =  1.0;
+    }
+
+    for(row=0;row<3;row++){
+      for(col=0;col<3;col++){
+        (*img)[i].rotMat[row][col] = 0.0;
+        for(j=0;j<3;j++)
+          (*img)[i].rotMat[row][col] += auxRotMat[row][j]*tempRotMat[j][col];
+      }
+    }
   }
 
-  /* Allocate moldata array */
+  /* Allocate moldata array.
+  */
   (*m)=malloc(sizeof(molData)*par->nSpecies);
   for( i=0; i<par->nSpecies; i++ ){
     (*m)[i].part = NULL;

src/lime.h

@@ -200,7 +200,7 @@ typedef struct {
   double freq,bandwidth;
   char *filename;
   double source_vel;
-  double theta,phi;
+  double theta,phi,incl,posang,azimuth;
   double distance;
   double rotMat[3][3];
 } image;

src/main.c

@@ -38,6 +38,7 @@ initParImg(inputPars *par, image **img)
   */
 
   int i,id,nImages;
+  const double defaultAngle=-999.0;
 
   /* Set 'impossible' default values for mandatory parameters */
   par->radius    = 0;
@@ -124,6 +125,9 @@ initParImg(inputPars *par, image **img)
     (*img)[i].molI=-1;
     (*img)[i].freq=-1.;
     (*img)[i].bandwidth=-1.;
+    (*img)[i].incl    = defaultAngle;
+    (*img)[i].azimuth = defaultAngle;
+    (*img)[i].posang  = defaultAngle;
   }
 
   /* Second-pass reading of the user-set parameters (this time just to read the par->moldatfile and img stuff). */