Added inclination/PA/azimuth to image parameters

[ghost]

Added the possibility to set the inclination, position angle and the azimuth angle in the image parameters (img[i].incl, img[i].posang, img[i].azimuth). If either of the three new angles are not set in the image structure the old theta and phi will be used.

[smaret]

For the record, I'm quoting @AttilaJuhasz 's email here:

I was a bit annoyed by the definition of the angles in the raytracer/image calculation so I changed that. Originally the theta and phi angles in the image definition denote the rotation angle around the x- and y-axes, which was a bit painful to calculate when modeling real observations. So I replaced theta and phi by the true inclination, position angle and azimuth angle, which are actually measured in real observations. Actually I left both version there so one can use either theta and phi or inclination, positiong angle and azimuth angle. The new angles are working fine unless one wants to have polarisation.

The issue/question on the polarisation is described in #68.

[smaret]

Are these really image parameters? It seems to me they describe the source, and not the way the images are computed. The only reason to put them in the img. structure would be to allow for different values in each image, but I don't think one would need to do that.

[ghost]

Inclination, position angle and azimuth are just as much image parameters as theta and phi. They are actually even more. Theta and phi describe the rotation around the X- and Y-axes and are equivalent to the inclination and position angles. However, it is more convenient to use inclination and position angle as observations traditionally report these angles not theta and phi. Moreover I don't see how you can describe the azimuth angle with theta and phi. For axisymmetric models the azimuth doesn't matter but for truly 3D models it does. Granted, you can always rotate the model structure with the azimuth but that's a bit troublesome as you need to implement the rotation in all model functions from the density() to the velocity(). It is way more simple to just rotate during the raytracing.

[smaret]

Inclination, position angle and azimuth are just as much image parameters as theta and phi.

Sorry I forgot that theta and phi were also defined in the img structure.

I've tested these changes with a 2D axis symmetric model and it seems to works fine. Indeed it is much more convenient to set the position and inclination angle in this fashion.

Could you please document these changes in the user manual? For this you need to edit the doc/usermanual.rst file. Fig. 1 (under doc/images) would probably need to be updated as well.

[smaret]

I would like to see this PR included in the next LIME release, but for this the new parameters needs to be documented.

[allegroLeiden]

This is a useful addition to LIME's capabilities. I found a couple of errors though, and I have some suggestions for enhancements.

• aux.c, line 194: where you have (*img)[i].phi do you mean (*img)[i].azimuth?

Also, having set defaults of -999.0, I would test all three angles for equality with that rather than testing for <-100. In this scheme, safest (read 'least unsafe') would be something like

const double defaultAngle=-999.0;
...
(*img)[i].azimuth=defaultAngle;
...etc
if(((*img)[i].azimuth==defaultAngle) || etc

Testing two floats for equality is evil, but the whole necessity for flagging unset optional parameters via impossible (or at least unlikely) default values is pretty evil, although unavoidable unless we want to force the user to set a bunch of extra boolean parameters. (Some day we will replace this C crap with a proper python interface which will make this stuff much neater!) At the very least I would test for <-900 rather than <-100, which is an angular value I could plausibly imagine a user choosing.

Finally, I would use || rather than | myself, even for expressions which will only take a non-zero value in the LSB. Bitwise operations should be reserved for bitwise contexts.

• aux.c, line 224: cos_incl here should be sin_incl.
• It seems to me that PA would be useful to have for the theta/phi setup as well (with default value 0). This might be a good stimulus to re-express the rotMat as a multiplication of individual matrices, which is much clearer and more robust. (I don't know now why I did not do that in the first place.) Suggested code (with some comment which I should also have included originally) would be:

    /*
The image rotation matrix is meant to rotate a point specified in an observer-anchored frame into the coordinate frame of the model. This is because it is primarily used in raytracing to calculate the position in the model frame of each ray starting point. In linear algebra terms, the model-frame position vector v_mod is related to the observer- (or image-) frame position vector v_obs 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_1 followed by R_2 followed by R_3, equation (2) can be expanded to give

    v_obs = R_3 * R_2 * R_1 * v_mod.            3

Inverting everything to return to the format of equation (1), which is what we need, we find

    v_mod = R_1^T * R_2^T * R_3^T * v_obs.      4

LIME provides two different schemes of {R_1, R_2, R_3}: {theta, phi, PA} and {azimuth, inclination, PA}. As an example, consider azimuth, which is a rotation of the model from the observer X axis towards the Y. This rotation matrix is thus

           ( cos(az) -sin(az) )
    R_az = (                  ).
           ( sin(az)  cos(az) )

However it is the transpose of this, R_az^T, which is used in the construction of R.
    */

    /* Load PA rotation matrix R_PA^T:
    */
    cos_pa   = cos((*img)[i].posang-PI/2.);
    sin_pa   = sin((*img)[i].posang-PI/2.);
    (*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;

    doThetaPhi = (((*img)[i].incl<-900.) || ((*img)[i].azimuth<-900.))?1:0;

    if(doThetaPhi){
        /* Load phi rotation matrix R_phi^T:
        */
        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^T:
        */
        cos_incl = cos((*img)[i].incl);
        sin_incl = sin((*img)[i].incl);
        auxRotMat[0][0] =  1.0;
        auxRotMat[0][1] =  0.0;
        auxRotMat[0][2] =  0.0;
        auxRotMat[1][0] =  0.0;
        auxRotMat[1][1] =  cos_incl;
        auxRotMat[1][2] =  sin_incl;
        auxRotMat[2][0] =  0.0;
        auxRotMat[2][1] = -sin_incl;
        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^T:
        */
        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^T:
        */
        cos_az   = cos((*img)[i].azimuth);
        sin_az   = sin((*img)[i].azimuth);
        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];
      }
    }

[allegroLeiden]

...or rather

           ( cos(az) -sin(az)  0 )
           (                     )
    R_az = ( sin(az)  cos(az)  0 ).
           (                     )
           (   0       0       1 )

[allegroLeiden]

After a great deal of offline discussion, I hope @AttilaJuhasz and I can now agree on the following:

(I adopt the usual convention that 0 represents an axis pointing out of the page and + represents one pointing into the page.)

1. That it is conservative programming practice to break the total rotation into components, one per axis.

2. That a rotation R of the model axes, seen from a fixed observer standpoint, has the same effect as a rotation R^T of the observer axes, seen from a fixed model standpoint (^T means transpose).

3. That if R = A*B*C then R^T = C^T * B^T * A^T, i.e. the components are applied in reverse order.

4. From the observer's point of view, the model XYZ axes transform under theta/phi rotations as follows:

{theta,phi} = {0,0}:

     ^ Y
     |
     |
<----+ Z
X

{theta,phi} = {90,0}:

<----+ Y
X    |
     |
     v Z

{theta,phi} = {0,90}:

     ^ Y
     |
     |
<----0 X
Z

{theta,phi} = {90,90}:

<----0 X
Y    |
     |
     v Z

5. That the correct order of {theta,phi} rotations, as seen from a fixed observer standpoint, is first theta, then phi, then PA.

6. That the construction of rotMat in Lime, which applies rotations of the observer axes from the fixed model standpoint, must therefore (from 2 and 3) invert this order, and apply the (appropriately transposed) rotation matrices in the order PA, phi, theta.

7. In other words, that the analysis of the present {theta,phi} setup, with the addition of PA, is exactly correct as I gave it on 1 July (with a caveat about PA which I will mention in a moment).

Moving on to the proposed {az,inc,PA} scheme, I believe that my analysis of @AttilaJuhasz's suggested rotMat (after the correction indicated) is correct. In particular:

8. The order of application of {az,inc,PA} from the fixed observer standpoint is az, then inc, then PA;

9. Thus from 2 and 3, when constructing rotMat this order must be inverted (=> as I have it in my 1 July).

I do suggest that all three of these angular values, as reverse-engineered from @AttilaJuhasz's supplied matrix, require either some applied scalar offset, or an inverse in sign, or both. This is purely a matter of convention and is something which astronomer users are in a better position to dictate than me; nevertheless, below is my understanding, subject to correction, of the way things ought to be arranged.

10. I would suggest that the desired orientation of model axes, seen from the fixed observer standpoint, when {az,inc,PA}={0,0,0}, is actually

   Z ^
     |
     |
   X 0---->
          Y

11. I would suggest that the direction of each of the three rotations from this point should be as diagrammed below:

{az,inc,PA} = {90,0,0}:

   Z ^
     |
     |
   Y +---->
          X

{az,inc,PA} = {0,90,0}:

   Z 0---->
     |    Y
     |
   X v

{az,inc,PA} = {0,0,90}:

   Y ^
     |
     |
<----0
Z    X

12. To get to the {0,0,0} starting point from the Lime definition of 'at rest' axes, we have to apply a 'pre-azimuth' of 90 and a 'pre-inclination' also of 90. Whichever values of az and inc the user specifies therefore, Lime should add 90deg to each before constructing the matrices.

13. Pseudo-code for az, inc and PA matrices, calculated from the observer-fixed standpoint, would therefore be as follows:

c = cos(az + PI/2.);
s = sin(az + PI/2.);

        (c -s  0)
azMat = (s  c  0);
        (0  0  1)

(cos(az+PI/2.) = -sin(az) etc, but let's keep it simple.)

c = cos(inc + PI/2.);
s = sin(inc + PI/2.);

         (1  0  0)
incMat = (0  c  s);
         (0 -s  c)

c = cos(PA);
s = sin(PA);

        ( c  s  0)
paMat = (-s  c  0);
        ( 0  0  1)

14. The az/inc/PA rotMat would thus be constructed as

rotMat = azMat^T * incMat^T * paMat^T

[ghost]

While I agree with most of what @allegroLeiden listed I do not agree with all of them.

4. For {theta, phi} = (90,0) the model Z-axis is pointed downwards in the image, anti-parallel with the declination axis. In most other modelling packages this is the other way around, i.e. an inclination of 90deg produces an image where the model Z-axis is pointed upwards. This should either be changed, or described clearly in the documentation as this can be confusing for full 3D models with no symmetry.

@allegroLeiden @smaret would it be to big of a change in the raytracer to move the observer at the starting position from z=-infinity to z=+infinity, so then for phi=0 and theta<90deg the model would face the observer with the side that is located at positive z values.

5)-6) The theta and phi rotations change both the inclination and the position angle in the general case. So I'm not really sure how to implement PA next to theta and phi in a self-consistent way, but maybe I just don't see the forest for the trees. :-) In the 1 July implementation of @allegroLeiden when using {theta, phi, PA} the first rotation by 'posang' will set the position angle of a disk in the X-Y plane of the model only if phi=0. For phi=90deg, the position angle of a disk in the x-y plane would be theta + posang if I'm not mistaken.

10. I think that for the {az, incl, PA} triplet, the most straightforward and cleanest implementation would probably be the following:

The observer would be located at z=+infinity looking down onto the model.
Starting orientation of the model axes as seen by the observer for {az, incl, PA} = (0, 0, 0) to be

Z 0---->
  |    Y
  |
X v

Then the three rotations would produce the following orientation

{az, incl, PA} = (90, 0, 0):

<----0 Z
X    |
     |
     v Y

{az, incl, PA} = (0, 90, 0):


Y +---->
  |    Z
  |
X v


{az, incl, PA} = (0, 0, 90):


Y ^
  |
  |
Z 0---->
       X

This way if we have a disk-like structure flattened along the z-axis, whose mid-plane corresponds to the model X-Y plane then {az,incl,PA}=(0,0,0) will give us an face-on orientation of the model, {az,incl,PA}=(0,90,0) would give us an image with an edge-on disk whose major axis is vertical on the image (i.e. aligned with the declination axis), while {az, incl, PA} = (0, 90, 90) would give us again an edge-on disk whose major axis is aligned horizontally in the image, which is consistent with the astronomical notation assuming that the RA axis of the image is pointed to the left.

[allegroLeiden]

• I have no opinion one way or the other about changes to the present theta/phi scheme, but if this is decided on, it should probably be put off until 2.0, because it would break backward compatibility. I'd encourage you to raise this as an issue; it should not be part of the present PR.
• PA is just a rotation about the observer's direction of gaze. As I understand the convention, positive values correspond to anticlockwise rotations. The sequence theta-phi-PA represents right-handed rotations successively about the observer X axis, observer Y axis then finally observer -ve Z axis. There is no rotation about the observer Z axis with just theta and phi alone.
• I am happy to accept any definition of az and incl but I think you ought to run your choice first not only by @smaret but also Michiel Hogerheijde.

[ghost]

1. Accepted and agreed.
2. You are right, PA is a rotation of the object around the line of sight. The definition used in accretion disks is that the PA of a projected disk, ring is the angle between the direction of North and the major axis of the disk, measured via East. Now lets take a thin ring in the X-Y planet centred on the origin of the model coordinate system. {theta, phi} = (>0deg, 0) will produce an image where the major axis of the projected ring will be aligned horizontally, i.e. parallel to the RA axis of the image. By definition this means PA=90deg. {theta, phi} = (0deg, >0deg) will also produce an inclined ring, but this time the major axis of the projected ring will be aligned vertically. This means that PA=0deg. So theta and phi do change the PA of the image.
3. Sure, absolutely.

[allegroLeiden]

2. Ok, I accept freely your definition of PA as it deals with astrophysical objects. In these terms, a combination of theta and phi can, I agree, give you PAs of multiples of pi/2. Two rotations alone will not, however, allow all possible orientations; a third is needed. Given that theta/phi is not, and perhaps was not intended to be, a system particularly well suited to astrophysical conventions, the simplest way to add the 3rd degree of freedom to me seems to be to add the Z rotation, which allows the additional simplicities as given in my code suggestion. It's irrelevant really whether one calls this Z rotation PA in the documentation or not. If there is not common agreement on this, though, and it seems there is not, then I would suggest leaving the theta/phi system just as it is, without adding the 3rd rotation, and just implementing your (very worthwhile) addition of az/incl/pa. Once you astronomers settle on the angle definitions, I am happy to work out the matrices. The one is not my field, the other most definitely is.

[smaret]

A figure in the user manual showing the various angle definition with respect to the x, y and z axis would be extremely useful I think.

[allegroLeiden]

There already is one for the theta/phi stuff.

[smaret]

Here is a figure showing the changes proposed by Ian, as I understand them.

Indeed the PA angle is usually defined from North to East, so we would need to change it sign. Also as noted by @AttilaJuhasz, an inclination of 0 degrees for disk rotating around the z axis usually corresponds to an edge-on disk; so the inclination should be replaced by 90-inclination in the figure above.

With these minor changes, the angle definitions would look like this:

Do we all agree on these definitions?

[allegroLeiden]

Sorry, @smaret , your diagrams don't seem to correspond with either @AttilaJuhasz 's original suggested code or his diagrams in his point 10 of 18 July. I have suggested no scheme myself, I am just trying to interpret AJ's. Without going into mathematical detail, your scheme and his cannot be consistent, because he assumes incl==0 corresponds to a face-on disk whereas you assume it corresponds to an edge-on disk.

[allegroLeiden]

I attach some drawings of a donut with its positive XY quadrant shaded, rotated to various angles according to AJ's scheme of 18 July (I'm assuming that his rotations are performed about the fixed frame and that they are done in the order az, incl, pa). @AttilaJuhasz please tell me if I have misunderstood. If not, the only question @smaret is whether you agree with this scheme or not.
az_in_pa007.pdf

[allegroLeiden]

Oops sorry it is a pdf.

[allegroLeiden]

Michiel agrees with your definitions @AttilaJuhasz. Pseudo code for matrices which will generate the behaviour you describe is

c = cos(az + PI/2.);
s = sin(az + PI/2.);

        ( c  s  0)
azMat = (-s  c  0);
        ( 0  0  1)

c = cos(inc + PI);
s = sin(inc + PI);

         (1  0  0)
incMat = (0  c  s);
         (0 -s  c)

c = cos(PA);
s = sin(PA);

        ( c  s  0)
paMat = (-s  c  0);
        ( 0  0  1)

rotMat = azMat^T * incMat^T * paMat^T

If you code the PR that way I'll accept it.

[smaret]

Superseded by #161.