mockgallib/src

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
#ifndef SKY_H
#define SKY_H 1

#include <vector>
#include <cstdio>
#include <cassert>
#include "halo.h"
#include "util.h"

class Sky {
 public:
  Sky(const double ra[], const double dec[], const double z[] );
  inline void compute_radec(const float x[], float* const radec) const;
  inline void compute_x(const float r, const float radec[], float* const x) const;
  double ra_range[2], dec_range[2], r_range[2];
  float left[3], right[3], width[3], centre[3]; // bounding box
  // r0, dec0 are the centre of ra, dec ranges, respectively.
  float ra0, dec0, theta0, cos_theta0, sin_theta0;
};

void Sky::compute_radec(const float x[], float* const radec) const
{
  // compute h->ra and dec from cuboid coordinate h->x
  // ra0, dec0: RA and Dec of xaxis (y=0, z=0)

  // theta measured from x-y plane
  float theta= asin(x[2]/util::norm(x));
    
  // rotate cosO in x-z plane
  float x1= cos_theta0*x[0] - sin_theta0*x[2];
  float y1= x[1];
  float rr= sqrt(x1*x1 + y1*y1);
  float phi= asin(y1/rr); if(x1 < 0) phi= M_PI - phi;

  radec[0]= ra0  - 180.0/M_PI*phi; // RA is right to left
  radec[1]= dec0 + 180.0/M_PI*theta;
    
  assert(-180.0f <= radec[1] && radec[1] <= 180.0f);
}

void Sky::compute_x(const float r, const float radec[], float* const x) const
{
  const float sinO = sin(radec[1]/180.0*M_PI);
  const float cosO = cos(radec[1]/180.0*M_PI);
  const float cos_phi = cos((ra0 - radec[0])/180.0*M_PI);
  const float sin_phi = sin((ra0 - radec[0])/180.0*M_PI);

  // usual mapping from r,ra,dec to x[3] with x-axis ra=ra0
  float x1[] = {r*cosO*cos_phi, r*cosO*sin_phi, r*sinO};

  // rotate x1[] theta0 about y axis
  x[0]=  cos_theta0*x1[0] + sin_theta0*x1[2];
  x[1]= x1[1];
  x[2]= -sin_theta0*x1[0] + cos_theta0*x1[2];
}


  
#endif

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
//
// Sky contains a volume information of ra/dec/r in the sky
//

#include <iostream>
#include <vector>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cassert>

#include "msg.h"
#include "cosmology.h"
#include "sky.h"


using namespace std;

static inline double rad(const double deg)
{
  return deg/180.0*M_PI;
}

static void compute_bounding_box(Sky * const sky);

Sky::Sky(const double ra[], const double dec[], const double z[])
{
  ra_range[0]= ra[0]; ra_range[1]= ra[1];
  dec_range[0]= dec[0]; dec_range[1]= dec[1];

  const double a_min= 1.0/(1.0 + z[0]);
  const double a_max= 1.0/(1.0 + z[1]);
  r_range[0]= cosmology_compute_comoving_distance(a_min);
  r_range[1]= cosmology_compute_comoving_distance(a_max);

  compute_bounding_box(this);

  // The x-axis (y=0, z=0) corresponds passes the centre of RA and Dec
  dec0= dec_range[0] + 0.5*(dec_range[1] - dec_range[0]);
  ra0=  ra_range[0] + 0.5*(ra_range[1] - ra_range[0]);

  theta0= rad(dec0);
  cos_theta0= cos(theta0);
  sin_theta0= sin(theta0);
}


static inline void set_cartisian(const double r, const double theta, const double phi, double * const x)
{
  // Return cartisian coordinate corresponds to polar coordinate
  // theta is measured from xy-plane, not from z axis
  const double rcosO= r*cos(theta);
  x[0]= rcosO*cos(phi);
  x[1]= rcosO*sin(phi);
  x[2]= r*sin(theta);
}

static inline void rotate_xz(double * const x, const double theta)
{
  // Rotate the vector by an angle theta (radian) in x-z plane
  const double cosO= cos(theta);
  const double sinO= sin(theta);
  
  double x1= cosO*x[0] - sinO*x[2];
  double z1= sinO*x[0] + sinO*x[2];

  x[0]= x1;
  x[2]= z1;
}

void compute_bounding_box(Sky* const sky)
{
  //
  // computes the size of cuboid that can enclose the sky region
  //
  // Output: sky->box[3]
  //
  const double theta= rad(0.5*fabs(sky->dec_range[1] - sky->dec_range[0]));
  const double phi= rad(0.5*(sky->ra_range[1] - sky->ra_range[0]));
  assert(phi> 0.0);
  assert(phi < 0.5*M_PI); // phi > pi/2 not implemented in this code

  const double theta_min= rad(sky->dec_range[0]);
  const double theta_max= rad(sky->dec_range[1]);
  const double theta0= 0.5*(theta_min + theta_max);
  double x[3];

  // x-range
  double x_min= sky->r_range[0];
  
  set_cartisian(sky->r_range[0], theta_min, phi, x);
  rotate_xz(x, -theta0);
  x_min= fmin(x_min, x[0]);
  
  set_cartisian(sky->r_range[0], theta_max, phi, x);
  rotate_xz(x, -theta0);
  x_min= fmin(x_min, x[0]);

  sky->right[0]= sky->r_range[1];
  sky->left[0]= x_min;
  sky->width[0]= sky->r_range[1] - x_min;

  // y-range
  double y_max= 0.0;
  
  set_cartisian(sky->r_range[1], theta_min, phi, x);
  rotate_xz(x, -theta0);
  y_max= fmax(y_max, x[1]);
  
  set_cartisian(sky->r_range[1], theta_max, phi, x);
  rotate_xz(x, -theta0);
  y_max= fmax(y_max, x[1]);

  sky->left[1]= -y_max;
  sky->right[1]= y_max;
  sky->width[1]= 2.0*y_max;

  // z-range
  sky->right[2]= sky->r_range[1]*sin(theta);
  sky->left[2]= -sky->right[2];
  sky->width[2]= 2.0*sky->right[2];

  sky->centre[0]= sky->left[0] + 0.5f*(sky->right[0] - sky->left[0]);
  sky->centre[1]= sky->centre[2]= 0.0f;
}