Skip to content

saclaymocks

Vendored geometry helpers and constants from SaclayMocks. Prefer installing SaclayMocks directly; this copy is a fallback.

box

box

Box sky<->cartesian geometry helpers (vendored from SaclayMocks).

Bundled copy of the SaclayMocks box module, used as a fallback when SaclayMocks is not installed. Provides the coordinate transforms between sky (ra, dec, R) and box-frame cartesian (X, Y, Z), plus box-geometry helpers. Prefer installing SaclayMocks directly.

ComputeXYZdeg 

ComputeXYZdeg(ra, dec, R, ra0, dec0)

Box cartesian coordinates from sky coordinates in degrees.

Degree-input wrapper around :func:ComputeXYZ.

Parameters:

Name Type Description Default
ra, dec array - like

Sky angles (degrees).

 required
R array - like

Comoving distance.

 required
ra0, dec0 float

Box-centre sky angles (degrees).

 required

Returns:

Name Type Description
tuple

(X, Y, Z) box-frame cartesian coordinates.
Source code in lelantos/saclaymocks/box.py 
15
16
17
18
19
20
21
22
23
24
25
26
27
28
def ComputeXYZdeg(ra,dec,R,ra0,dec0) :
    """Box cartesian coordinates from sky coordinates in degrees.

    Degree-input wrapper around :func:`ComputeXYZ`.

    Args:
        ra, dec (array-like): Sky angles (degrees).
        R (array-like): Comoving distance.
        ra0, dec0 (float): Box-centre sky angles (degrees).

    Returns:
        tuple: ``(X, Y, Z)`` box-frame cartesian coordinates.
    """
    return ComputeXYZ(np.radians(ra),np.radians(dec),R,np.radians(ra0),np.radians(dec0))

ComputeXYZ 

ComputeXYZ(ra, dec, R, ra0, dec0)

XYZ of a point P (ra,dec,R) in a frame with observer at O, Z along OP, X along ra0, Y along dec0

angles in radians tested that ra,dec, R = box.ComputeRaDecR(R0,ra0,dec0,X,Y,Z) x,y,z = box.ComputeXYZ(ra[0],dec[0],R,ra0,dec0) print(x-X,y-Y,z-R0-Z)# prints ~1E-13 for random inputs

Source code in lelantos/saclaymocks/box.py 
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
def ComputeXYZ(ra,dec,R,ra0,dec0) :
    '''
    XYZ of a point P (ra,dec,R) in a frame with
    observer at O, Z along OP, X along ra0, Y along dec0

    angles in radians
    tested that ra,dec, R = box.ComputeRaDecR(R0,ra0,dec0,X,Y,Z)
    x,y,z = box.ComputeXYZ(ra[0],dec[0],R,ra0,dec0)
    print(x-X,y-Y,z-R0-Z)#        prints ~1E-13  for random inputs
    '''

    theta0 = PI/2 - dec0 # polar angle or colatitude
                #......         unit vector in X, Y and Z directions
    cosra0 = np.cos(ra0)
    sinra0 = np.sin(ra0)
    costheta0 = np.cos(theta0)
    sintheta0 = np.sin(theta0)
    xVec = np.array([-sinra0, cosra0 ,0])
    yVec = np.array([-costheta0 * cosra0, -costheta0 * sinra0 , sintheta0])
    zVec = np.array([sintheta0 * cosra0, sintheta0 * sinra0, costheta0])
                # coordinates of OP
    theta = PI/2 - dec
    cosra = np.cos(ra)
    sinra = np.sin(ra)
    costheta = np.cos(theta)
    sintheta = np.sin(theta)
    OP = R * np.array([sintheta*cosra,sintheta*sinra,costheta])

    #print('op',OP)
    #print(xVec,yVec,zVec)
    X = np.dot(OP,xVec)
    Y = np.dot(OP,yVec)
    Z = np.dot(OP,zVec)
    return X,Y,Z

ComputeRaDecR 

ComputeRaDecR(R0, ra0, dec0, X, Y, Z)

ra, dec, R of P(X,Y,Z) in a box with center located in C(R0,ra0,dec0)

angles in radian -pi < ra < pi -pi/2 < dec < pi/2 R0, ra0, dec0 define the position of the box in a (U,V,W) frame centered at the oberver position O box oriented such that Z along OC, X along ra, Y along dec R0,ra0,dec0 are scalar X, Y, Z can be scalars, 1D or 2D arrays

Source code in lelantos/saclaymocks/box.py 
 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
def ComputeRaDecR(R0,ra0,dec0,X,Y,Z) :
    '''
    ra, dec, R of P(X,Y,Z) in a box with center located in C(R0,ra0,dec0)

    angles in radian -pi < ra < pi   -pi/2 < dec < pi/2
    R0, ra0, dec0 define the position of the box in a (U,V,W) frame
        centered at the oberver position O
    box oriented such that Z along OC, X along ra, Y along dec
    R0,ra0,dec0 are scalar   X, Y, Z can be scalars, 1D or 2D arrays
    '''
    theta0 = np.pi/2 - dec0 # polar angle or colatitude
                #......         unit vector in X, Y and Z directions
    cosra0 = np.cos(ra0)
    sinra0 = np.sin(ra0)
    costheta0 = np.cos(theta0)
    sintheta0 = np.sin(theta0)
    # unit vector along X Y and Z in the (U,V,W) frame
    xVec = np.array([-sinra0, cosra0 ,0])[:, None,None] # (3,1,1)
    yVec = np.array([-costheta0 * cosra0, -costheta0 * sinra0 , sintheta0])[:, None,None]
    zVec = np.array([sintheta0 * cosra0, sintheta0 * sinra0, costheta0])[:, None,None]
    #print(xVec.shape, yVec.shape, zVec.shape)
    #print(X.shape, Y.shape, Z.shape)
                #......         vector OP
    OPVec = X * xVec + Y * yVec + (R0+Z) * zVec # broadcasting due to [:, None,None] above
                    # allows to run with X and Y that are 1D or 2D arrays
    U = OPVec[0]    # (1,1) (1,N) or (N1,N2) depending on X,Y,Z shape
    V = OPVec[1]
    W = OPVec[2]
                #......         ra dec of OP
                # arctan2(y,x): angle relative to x axis in [-pi,pi]
    ra = np.arctan2(V,U)    # angle raltive to U
    dec = np.arctan2(W, np.sqrt(U*U + V*V))  # angle relative to (u,V) plane
                                        # in [-pi2/,pi/2] since sqrt(U*U+V*V) > 0
    R = np.sqrt(U*U+V*V+W*W)
    #print(ra.shape[0],ra.shape[1])
    if (np.isscalar(X+Y+Z)):
        return ra[0,0],dec[0,0],R[0,0]
    #if (np.isscalar((X+Y+Z)[0])) : # X,Y,Z 1D arrays
    #   we may get X+Y+Z a zero length array, in which case, this fails
    if (ra.shape[0]==1) : # X,Y,Z 1D arrays
        return ra[0],dec[0],R[0]
    return ra, dec, R

Compute_cos_min 

Compute_cos_min(LX, LY, Rmax)

Maximum tangents and minimum cosine of the box opening cone.

Parameters:

Name Type Description Default
LX, LY float

Transverse box sizes.

 required
Rmax float

Maximum comoving distance.

 required

Returns:

Name Type Description
tuple

(tanx_max, tany_max, cos_min).
Source code in lelantos/saclaymocks/box.py 
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
def Compute_cos_min(LX,LY,Rmax) :
    """Maximum tangents and minimum cosine of the box opening cone.

    Args:
        LX, LY (float): Transverse box sizes.
        Rmax (float): Maximum comoving distance.

    Returns:
        tuple: ``(tanx_max, tany_max, cos_min)``.
    """
    #
    #  angle for a point at X=Xmax and Rmax from observer
    sinx_max = (LX/2) / Rmax
    cosx_min = np.sqrt(1 - sinx_max**2)
    tanx_max = sinx_max / cosx_min
    #  angle for a point at Y=Ymax and Rmax from observer
    siny_max = (LY/2) / Rmax
    cosy_min = np.sqrt(1 - siny_max**2)
    tany_max = siny_max / cosy_min
    #
    sin_max = (np.sqrt(LX*LX+LY*LY)/2) / Rmax
    cos_min = np.sqrt(1 - sin_max**2)
    return tanx_max,tany_max,cos_min

box_center 

box_center(LX, LY, Rmax, Rmin, margin)

Radial box centre and depth enclosing a cone within a transverse box.

Parameters:

Name Type Description Default
LX, LY float

Transverse box sizes.

 required
Rmax, Rmin float

Comoving distance bounds of the cone.

 required
margin float

Margin kept inside the box edges.

 required

Returns:

Name Type Description
tuple

(R0, LZ, tanx_max, tany_max) — box-centre distance, depth and

opening tangents.
Source code in lelantos/saclaymocks/box.py 
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
def box_center(LX,LY,Rmax,Rmin,margin) :
    """Radial box centre and depth enclosing a cone within a transverse box.

    Args:
        LX, LY (float): Transverse box sizes.
        Rmax, Rmin (float): Comoving distance bounds of the cone.
        margin (float): Margin kept inside the box edges.

    Returns:
        tuple: ``(R0, LZ, tanx_max, tany_max)`` — box-centre distance, depth and
        opening tangents.
    """
    # we want to include a cone Rmin < R < Rmax in a LX x LY box minus margin
    # returns the resulting tan_max, the required LZ and the box_center distance
    tanx_max,tany_max,cos_min = Compute_cos_min(LX-2*margin,LY-2*margin,Rmax)
    LZ = Rmax - (cos_min * Rmin) + 2 * margin
    R0 = Rmax -LZ/2 + margin
    return R0,LZ,tanx_max,tany_max

box_limitOld 

box_limitOld(LX, LY, LZ, R0, margin)

Legacy: radial/angular limits of a box at distance R0 (see box_limit).

Parameters:

Name Type Description Default
LX, LY, LZ float

Box sizes.

 required
R0 float

Box-centre comoving distance.

 required
margin float

Margin kept inside the box edges.

 required

Returns:

Name Type Description
tuple

(Rmin, Rmax, tanx_max, tany_max).
Source code in lelantos/saclaymocks/box.py 
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
def box_limitOld(LX,LY,LZ,R0,margin) :
    """Legacy: radial/angular limits of a box at distance R0 (see box_limit).

    Args:
        LX, LY, LZ (float): Box sizes.
        R0 (float): Box-centre comoving distance.
        margin (float): Margin kept inside the box edges.

    Returns:
        tuple: ``(Rmin, Rmax, tanx_max, tany_max)``.
    """
    # compute R_min, R_max, tanx_max and tany_max
    # to be inside the box with a margin
    # for a box located at R0 from the observer
    #    we want that in any direction R_min < R_QSO < R_max
    Rmax = R0 + LZ/2  - margin
    #  angle for a point at X=Xmax and Rmax from observer
    sinx_max = (LX/2 -margin) / Rmax
    cosx_min = np.sqrt(1 - sinx_max**2)
    tanx_max = sinx_max / cosx_min
#    print(sinx_max, tanx_max)    # prov
    siny_max = (LY/2 -margin) / Rmax
    cosy_min = np.sqrt(1 - siny_max**2)
    tany_max = siny_max / cosy_min
    cos_min = min(cosx_min,cosy_min)
#    sin_max = max(sinx_max,siny_max)
#    cos_min = np.sqrt(1 - sin_max**2)
    Rmin = (R0 - LZ/2 + margin) / cos_min

#    print(R0-LZ/2, Rmin, Rmax)
    #print(0, Rmin-R0+LZ/2, Rmax-R0+LZ/2)
    #print(sin_alpha)
    return Rmin,Rmax,tanx_max,tany_max

box_limit 

box_limit(LX, LY, LZ, R0, margin)

Radial and angular limits of a box at distance R0 (with a margin).

Parameters:

Name Type Description Default
LX, LY, LZ float

Box sizes (requires LX == LY).

 required
R0 float

Box-centre comoving distance.

 required
margin float

Margin kept inside the box edges.

 required

Returns:

Name Type Description
tuple

(Rmin, Rmax, tanx_max, tany_max).
Source code in lelantos/saclaymocks/box.py 
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
def box_limit(LX,LY,LZ,R0,margin) :
    """Radial and angular limits of a box at distance R0 (with a margin).

    Args:
        LX, LY, LZ (float): Box sizes (requires ``LX == LY``).
        R0 (float): Box-centre comoving distance.
        margin (float): Margin kept inside the box edges.

    Returns:
        tuple: ``(Rmin, Rmax, tanx_max, tany_max)``.
    """
    # compute R_min, R_max, tanx_max and tany_max
    # to be inside the box with a margin
    # for a box located at R0 from the observer
    #    we want that in any direction R_min < R_QSO < R_max
    Rmax = R0 + LZ/2  - margin
    #  angle for a point at X=Xmax and Rmax from observer
    if (LX != LY) :
        print("case LX =",LX,"!= LY=",LY,"not implemented yet")
        exit(0)
    sinx_max = (LX/2 -margin) / Rmax
    cosx_min = np.sqrt(1 - sinx_max**2)
    tanx_max = sinx_max / cosx_min
#    print(sinx_max, tanx_max)    # prov
    siny_max = (LY/2 -margin) / Rmax
    cosy_min = np.sqrt(1 - siny_max**2)
    tany_max = siny_max / cosy_min

    sin_max = (np.sqrt(LX*LX+LY*LY)/2 -margin) / Rmax
    cos_min = np.sqrt(1 - sin_max**2)
    Rmin = (R0 - LZ/2 + margin) / cos_min

#    print(R0-LZ/2, Rmin, Rmax)
    #print(0, Rmin-R0+LZ/2, Rmax-R0+LZ/2)
    #print(sin_alpha)
    return Rmin,Rmax,tanx_max,tany_max

sample_box 

sample_box(xbox, ybox, zbox, rho, threshold=2.5)

Obsolete: sample QSO positions at density peaks of a box.

Parameters:

Name Type Description Default
xbox, ybox, zbox array - like

Box cell-centre coordinates per axis.

 required
rho ndarray

Density field.

 required
threshold float

Peak threshold in units of the rms.

 2.5

Returns:

Name Type Description
list

The sampled QSO objects.
Source code in lelantos/saclaymocks/box.py 
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
def sample_box(xbox,ybox,zbox,rho,threshold=2.5):
    """Obsolete: sample QSO positions at density peaks of a box.

    Args:
        xbox, ybox, zbox (array-like): Box cell-centre coordinates per axis.
        rho (numpy.ndarray): Density field.
        threshold (float, optional): Peak threshold in units of the rms.

    Returns:
        list: The sampled QSO objects.
    """
#                                                               obsolete
#
    dmax = 3    # we compute rho over +- dmax cells  should be a parameter <==

#.......................    select pixels with rho > threshold*rms
    rms = rho.std()
    iqso = sp.where(rho > threshold*rms)
#    print("threshold*rms=", threshold*rms)
    iqso = sp.array(iqso)        # (3,N_QSO) 3 raws N col
#.......................   remove QSO at less than dmax cells from the box edges
    imax = xbox.size - dmax
    jmax = ybox.size - dmax
    kmax = zbox.size - dmax
    jqso = iqso.T
    jqso = jqso[jqso[:,0]>dmax]
    jqso = jqso[jqso[:,0]<imax]
    jqso = jqso[jqso[:,1]>dmax]
    jqso = jqso[jqso[:,1]<jmax]
    jqso = jqso[jqso[:,2]>dmax]
    jqso = jqso[jqso[:,2]<kmax]
    iqso = jqso.T
#    print(iqso)
#    print(sp.where(iqso))

#   needs to implement z distribution
#   taking into account that z depends not only on Z but also on X and Y
#   must also select in a cone, not in the full box, this requires to know z0
#   must randomize the position inside the cell

    print("found: ",len(iqso[0])," quasars")

#    return [x+y+z for x,y,z in zip(xbox[iqso[0]],ybox[iqso[1]],zbox[iqso[2]])]
    return [qso(x,y,z) for x,y,z in zip(xbox[iqso[0]],ybox[iqso[1]],zbox[iqso[2]])]

ComputeRaDecR2 

ComputeRaDecR2(x, y, z, ra0, dec0)

Sky coordinates from box cartesian coordinates (rotation form).

Parameters:

Name Type Description Default
x, y, z array - like

Box-frame cartesian coordinates.

 required
ra0, dec0 float

Box-centre sky angles (radians).

 required

Returns:

Name Type Description
tuple

(ra, dec, R) with angles in radians.
Source code in lelantos/saclaymocks/box.py 
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
def ComputeRaDecR2(x, y, z, ra0, dec0):
    """Sky coordinates from box cartesian coordinates (rotation form).

    Args:
        x, y, z (array-like): Box-frame cartesian coordinates.
        ra0, dec0 (float): Box-centre sky angles (radians).

    Returns:
        tuple: ``(ra, dec, R)`` with angles in radians.
    """
    # all angles should be in radians
    # x, y, z are arrays of same dimension, and ra0 dec0 floats
    # Compute (ra,dec,R0) in base R of a point P(x,y,z) where (x,y,z)
    # are the coordinates in a base R' rotated by ra0 and dec0 with
    # respect to base R :
    # P_in_R' = Rot(dec0).Rot(ra0).P_in_R
    numra = (np.cos(ra0)*x - np.sin(dec0)*np.sin(ra0)*y
             + np.cos(dec0)*np.sin(ra0)*z)
    denomra = (-np.sin(ra0)*x - np.sin(dec0)*np.cos(ra0)*y
               + np.cos(dec0)*np.cos(ra0)*z)
    numdec = np.cos(dec0)*y + np.sin(dec0)*z
    R = np.sqrt(x**2 + y**2 + z**2)

    ra = np.zeros(x.shape)
    msk = np.where((numra > 0) & (denomra > 0))
    ra[msk] = np.arctan(numra[msk] / denomra[msk])
    msk = np.where((numra > 0) & (denomra < 0))
    ra[msk] = np.arctan(numra[msk] / denomra[msk]) + np.pi
    msk = np.where((numra < 0) & (denomra < 0))
    ra[msk] = np.arctan(numra[msk] / denomra[msk]) + np.pi
    msk = np.where((numra < 0) & (denomra > 0))
    ra[msk] = np.arctan(numra[msk] / denomra[msk]) + 2*np.pi

    msk = np.where((numra == 0) & (denomra > 0))
    ra[msk] = 0
    msk = np.where((numra > 0) & (denomra == 0))
    ra[msk] = np.pi / 2
    msk = np.where((numra == 0) & (denomra < 0))
    ra[msk] = np.pi
    msk = np.where((numra < 0) & (denomra == 0))
    ra[msk] = 3*np.pi / 2
    msk = np.where((numra == 0) & (denomra == 0))
    ra[msk] = 0

    dec = np.arcsin(numdec/R)
    return ra, dec, R

ComputeXYZ2 

ComputeXYZ2(ra, dec, R, ra0, dec0)

Box cartesian coordinates from sky coordinates (rotation form).

Parameters:

Name Type Description Default
ra, dec array - like

Sky angles (radians).

 required
R array - like

Comoving distance.

 required
ra0, dec0 float

Box-centre sky angles (radians).

 required

Returns:

Name Type Description
tuple

(x, y, z) box-frame cartesian coordinates.
Source code in lelantos/saclaymocks/box.py 
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
def ComputeXYZ2(ra, dec, R, ra0, dec0):
    """Box cartesian coordinates from sky coordinates (rotation form).

    Args:
        ra, dec (array-like): Sky angles (radians).
        R (array-like): Comoving distance.
        ra0, dec0 (float): Box-centre sky angles (radians).

    Returns:
        tuple: ``(x, y, z)`` box-frame cartesian coordinates.
    """
    # all angles should be in radians
    x = R * (np.cos(ra0)*np.cos(dec)*np.sin(ra)
             - np.sin(ra0)*np.cos(dec)*np.cos(ra))
    y = R * (-np.sin(ra0)*np.sin(dec0)*np.cos(dec)*np.sin(ra)
             + np.cos(dec0)*np.sin(dec)
             - np.cos(ra0)*np.sin(dec0)*np.cos(dec)*np.cos(ra))
    z = R * (np.cos(dec0)*np.sin(ra0)*np.cos(dec)*np.sin(ra)
             + np.sin(dec0)*np.sin(dec)
             + np.cos(ra0)*np.cos(dec0)*np.cos(dec)*np.cos(ra))
    return x, y, z

constant

constant

Physical constants and fiducial parameters (vendored from SaclayMocks).

Bundled copy of the SaclayMocks constant module, used as a fallback when SaclayMocks is not installed. Prefer installing SaclayMocks directly.