Skip to content

Commit 48f00a9

Browse files
slizewskslizewsk
committed
Adding main functions
1 parent 5224504 commit 48f00a9

File tree

15 files changed

+1896
-20
lines changed

15 files changed

+1896
-20
lines changed

‎.DS_Store‎

8 KB
Binary file not shown.

‎.github/.DS_Store‎

6 KB
Binary file not shown.

‎.github/workflows/ci.yml‎

Lines changed: 30 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,30 @@
1+
name: Run Tests
2+
3+
on:
4+
push:
5+
branches:
6+
- main
7+
pull_request:
8+
9+
jobs:
10+
test:
11+
runs-on: ubuntu-latest
12+
13+
steps:
14+
- name: Checkout repo
15+
uses: actions/checkout@v4
16+
17+
- name: Set up Python
18+
uses: actions/setup-python@v4
19+
with:
20+
python-version: "3.11"
21+
22+
- name: Install dependencies
23+
run: |
24+
python -m pip install --upgrade pip
25+
pip install pytest
26+
pip install -e '.[test]'
27+
28+
- name: Run tests
29+
run: pytest
30+
# run: pytest df_sampling/tests/test_sampling.py

‎.github/workflows/sub_package_update.yml‎

Lines changed: 7 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -8,13 +8,13 @@ permissions:
88
on:
99
# Allow manual runs through the web UI
1010
workflow_dispatch:
11-
schedule:
12-
# ┌───────── minute (0 - 59)
13-
# │ ┌───────── hour (0 - 23)
14-
# │ │ ┌───────── day of the month (1 - 31)
15-
# │ │ │ ┌───────── month (1 - 12 or JAN-DEC)
16-
# │ │ │ │ ┌───────── day of the week (0 - 6 or SUN-SAT)
17-
- cron: '0 7 * * 1' # Every Monday at 7am UTC
11+
# schedule:
12+
# # ┌───────── minute (0 - 59)
13+
# # │ ┌───────── hour (0 - 23)
14+
# # │ │ ┌───────── day of the month (1 - 31)
15+
# # │ │ │ ┌───────── month (1 - 12 or JAN-DEC)
16+
# # │ │ │ │ ┌───────── day of the week (0 - 6 or SUN-SAT)
17+
# - cron: '0 7 * * 1' # Every Monday at 7am UTC
1818

1919
jobs:
2020
update:

‎df_sampling/.DS_Store‎

6 KB
Binary file not shown.

‎df_sampling/__init__.py‎

Lines changed: 9 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,11 @@
1-
from .example_mod import do_primes
1+
from .dosampling import Params, ParamsHernquist, DataSampler
2+
from .make_obs import mockobs
3+
24
from .version import version as __version__
35

4-
# Then you can be explicit to control what ends up in the namespace,
5-
__all__ = ['do_primes']
6+
__all__ = [
7+
"Params",
8+
"ParamsHernquist",
9+
"DataSampler",
10+
"mockobs",
11+
]

‎df_sampling/coord_transforms.py‎

Lines changed: 283 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,283 @@
1+
from .core_imports import np
2+
3+
def cart2sph_pos(xyz):
4+
x,y,z = xyz
5+
D = np.sqrt(x**2+y**2) # projected dist
6+
r = np.sqrt(x**2 + y**2 + z**2) #radial dist
7+
phi = np.arctan2(y, x) #azimuthal angle
8+
theta = np.arctan(z / D) #polar angle
9+
# if phi < 0:
10+
# phi += 2 * np.pi
11+
# same as
12+
# theta = np.arcsin(z / r) #polar angle
13+
return np.array([r,theta,phi])
14+
def sph2cart_pos(rtp):
15+
r,theta,phi = rtp
16+
x = r * np.cos(theta) * np.cos(phi)
17+
y = r * np.cos(theta) * np.sin(phi)
18+
z = r * np.sin(theta)
19+
return np.array([x, y, z])
20+
21+
def cart2sph_vel(vel_cart,xyz):
22+
x,y,z = xyz
23+
vx,vy,vz = vel_cart
24+
dist,theta,phi = cart2sph_pos(xyz);
25+
proj_dist = np.sqrt(x**2 + y**2)
26+
vr = np.dot(xyz,vel_cart)/dist
27+
mu_theta = ((z * (x * vx + y * vy) - np.square(proj_dist) * vz)) / np.square(dist) / proj_dist
28+
vtheta = -mu_theta * dist
29+
30+
mu_phi = (x * vy - y * vx) / np.square(proj_dist)
31+
vphi = mu_phi * dist * np.cos(theta)
32+
33+
vel_sph = np.array([vr,vtheta,vphi])
34+
return vel_sph
35+
def sph2cart_vel(vel_sph,rtp):
36+
r, theta,phi = rtp
37+
sph_mat = np.array([[np.cos(phi) * np.cos(theta), -np.cos(phi) * np.sin(theta), -np.sin(phi)],
38+
[np.sin(phi) * np.cos(theta), -np.sin(phi) * np.sin(theta), np.cos(phi)],
39+
[np.sin(theta), np.cos(theta), 0] ])
40+
vel_cart = np.dot(sph_mat,vel_sph)
41+
return vel_cart
42+
43+
def eq2gc_pos(r_eq,
44+
R =np.array([[-0.05487395617553902, -0.8734371822248346,-0.48383503143198114],
45+
[0.4941107627040048, -0.4448286178025452,0.7469819642829028],
46+
[-0.8676654903323697, -0.1980782408317943,0.4559842183620723]]) ,
47+
H = np.array([[0.9999967207734917, 0.0, 0.002560945579906427],
48+
[0.0, 1.0, 0.0],
49+
[-0.002560945579906427, 0.0, 0.9999967207734917]]),
50+
offsett = np.array([8.112,0,0.0208]),return_cart=False):
51+
dist,ra,dec = r_eq
52+
# print('pos_eq',r_eq)
53+
r_icrs = sph2cart_pos(np.array([dist,dec,ra]))
54+
# print('xyz_helio',r_icrs)
55+
r_gal = np.dot(R,r_icrs)
56+
# print('rot1',r_gal)
57+
r_gal -= offsett
58+
r_gal = np.dot(H ,r_gal)
59+
if return_cart:
60+
return r_gal
61+
# print('gc xyz',r_gal)
62+
else:
63+
r,theta,phi = cart2sph_pos(r_gal)
64+
return np.array([r,theta,phi])
65+
def gc2eq_pos(xyz,
66+
R = np.array([[-0.05487395617553902, -0.8734371822248346,-0.48383503143198114],
67+
[0.4941107627040048, -0.4448286178025452,0.7469819642829028],
68+
[-0.8676654903323697, -0.1980782408317943,0.4559842183620723]]) ,
69+
H = np.array([[0.9999967207734917, 0.0, 0.002560945579906427],
70+
[0.0, 1.0, 0.0],
71+
[-0.002560945579906427, 0.0, 0.9999967207734917]]),
72+
offsett = np.array([8.112,0,0.0208]),return_cart=False):
73+
# print('gc xyz',xyz)
74+
H_inv = np.linalg.inv(H)
75+
R_inv = np.linalg.inv(R)
76+
# H_inv = np.dot(np.array([[-1, 0, 0],[0, 1, 0],[0, 0, -1]]),H)
77+
# R_inv = np.dot(np.array([[-1, 0, 0],[0, 1, 0],[0, 0, -1]]),R)
78+
r_gal = np.dot(H_inv,xyz)
79+
r_gal += offsett
80+
# print('rot1',r_gal)
81+
r_icrs = np.dot(R_inv, r_gal)
82+
if return_cart:
83+
return r_icrs
84+
else:
85+
r,dec,ra = cart2sph_pos(r_icrs)
86+
return np.array([r,ra,dec])
87+
88+
def eq2gc_vel(vel_eq,pos_eq,
89+
R =np.array([[-0.05487395617553902, -0.8734371822248346,-0.48383503143198114],
90+
[0.4941107627040048, -0.4448286178025452,0.7469819642829028],
91+
[-0.8676654903323697, -0.1980782408317943,0.4559842183620723]]) ,
92+
H = np.array([[0.9999967207734917, 0.0, 0.002560945579906427],
93+
[0.0, 1.0, 0.0],
94+
[-0.002560945579906427, 0.0, 0.9999967207734917]]),
95+
offsett = np.array([8.112,0,0.0208]),
96+
solarmotion=np.array([12.9/100, 245.6/100, 7.78/100]),
97+
return_cart=False):
98+
dist,ra,dec = pos_eq
99+
vlos,pmra,pmdec = vel_eq
100+
vlos /= 100
101+
conversion_factor = 4.740470463533349
102+
dist_with_conversion = dist * conversion_factor
103+
104+
vra = (dist_with_conversion * pmra ) / 100
105+
vdec = (dist_with_conversion * pmdec) / 100
106+
107+
r_gal = eq2gc_pos(pos_eq,return_cart=True)
108+
v_icrs = sph2cart_vel(rtp=np.array([dist,dec,ra ]), vel_sph=np.array([vlos, vdec, vra]))
109+
A = H @ R
110+
v_gal = A @ v_icrs
111+
v_gal += solarmotion
112+
v_gal_sph = cart2sph_vel(xyz=r_gal,vel_cart=v_gal)
113+
if return_cart:
114+
return v_gal
115+
else:
116+
return v_gal_sph
117+
def gc2eq_vel(vel_gc,pos_gc,
118+
R =np.array([[-0.05487395617553902, -0.8734371822248346,-0.48383503143198114],
119+
[0.4941107627040048, -0.4448286178025452,0.7469819642829028],
120+
[-0.8676654903323697, -0.1980782408317943,0.4559842183620723]]) ,
121+
H = np.array([[0.9999967207734917, 0.0, 0.002560945579906427],
122+
[0.0, 1.0, 0.0],
123+
[-0.002560945579906427, 0.0, 0.9999967207734917]]),
124+
offsett = np.array([8.112,0,0.0208]),solarmotion= np.array([12.9/100, 245.6/100, 7.78/100]),
125+
in_cart=False,return_cart=False):
126+
if in_cart:
127+
ricrs = pos_gc
128+
v_gal = vel_gc
129+
else:
130+
ricrs = sph2cart_pos(pos_gc)
131+
v_gal = sph2cart_vel(vel_sph=vel_gc,rtp=pos_gc)
132+
#
133+
dist,ra,dec = gc2eq_pos(ricrs, return_cart=False)
134+
xyz_helio = sph2cart_pos(np.array([dist,dec,ra]))
135+
# print('distance',dist)
136+
H_inv = np.linalg.inv(H)
137+
R_inv = np.linalg.inv(R)
138+
# input spherical velocity and position
139+
# print('v_gal',v_gal)
140+
v_gal -= solarmotion
141+
A = np.dot(R_inv,H_inv)
142+
# v_icrs = np.linalg.inv(A) @ v_gal
143+
v_icrs = A @ v_gal
144+
# print('v_icrs',v_icrs)
145+
if return_cart:
146+
return v_icrs*100
147+
else:
148+
vlos,vtheta,vphi = cart2sph_vel(vel_cart = v_icrs, xyz = xyz_helio)
149+
# vlos,vtheta,vphi = cart2sph_vel(vel_cart = np.array([v_icrs[0],v_icrs[1],v_icrs[2]]),
150+
# xyz = np.array([xyz_helio[0],xyz_helio[1],xyz_helio[2]]))
151+
152+
conversion_factor = 4.740470463533349
153+
dist_with_conversion = dist * conversion_factor
154+
# print('vra',vphi)
155+
vtheta *= 100
156+
vphi *= 100
157+
vlos *= 100
158+
pmra = (vphi) / dist_with_conversion
159+
pmdec = (vtheta) / dist_with_conversion
160+
vel_eq = np.array([vlos,pmra,pmdec])
161+
# print("v_eq:", vel_eq)
162+
return vel_eq
163+
164+
165+
166+
# # defining transformation matrices to ref/call
167+
# R = np.array([[-0.05487395617553902, -0.8734371822248346, -0.48383503143198114],
168+
# [0.4941107627040048, -0.4448286178025452, 0.7469819642829028],
169+
# [-0.8676654903323697, -0.1980782408317943, 0.4559842183620723]])
170+
# H = np.array([[0.9999967207734917, 0.0, 0.002560945579906427],
171+
# [0.0, 1.0, 0.0],
172+
# [-0.002560945579906427, 0.0, 0.9999967207734917]])
173+
# H_inv = np.dot(np.array([[1, 0, 0],[0, 1, 0],[0, 0, 1]]),H)
174+
# R_inv = np.dot(np.array([[1, 0, 0],[0, 1, 0],[0, 0, 1]]),R)
175+
176+
# these four are from jeffs code
177+
# transforms cartesian <-> spherical positions or velocities
178+
# def c2s_pos(x, y, z):
179+
# dist = np.sqrt(np.dot(np.array([x, y, z]),np.array([x, y, z])))
180+
# phi = np.arctan2(y, x)
181+
# theta = np.arcsin(z / dist)
182+
# return np.array([dist, theta, phi])
183+
# def s2c_pos(r, theta, phi):
184+
# x = r * np.cos(phi) * np.cos(theta)
185+
# y = r * np.sin(phi) * np.cos(theta)
186+
# z = r * np.sin(theta)
187+
# return np.array([x, y, z])
188+
189+
# def c2s_vel(x, y, z, vx, vy, vz):
190+
# sph_pos = c2s_pos(x, y, z);
191+
# dist = sph_pos[0]
192+
# lat = sph_pos[1]
193+
# lon = sph_pos[2]
194+
195+
# proj_dist = np.sqrt(np.square(x) + np.square(y))
196+
197+
# vr = np.dot(np.array([x, y, z]), np.array([vx, vy, vz])) / dist
198+
199+
# mu_theta = ((z * (x * vx + y * vy) -
200+
# np.square(proj_dist) * vz)
201+
# ) / np.square(dist) / proj_dist
202+
# vtheta = -mu_theta * dist
203+
204+
# mu_phi = (x * vy - y * vx) / np.square(proj_dist)
205+
# vphi = mu_phi * dist * np.cos(lat)
206+
207+
# return np.array([vr, vtheta, vphi])
208+
# def s2c_vel(r, theta, phi, vr, vtheta,vphi):
209+
# vx = vr * np.cos(phi) * np.cos(theta) - vphi * np.sin(phi)- vtheta * np.cos(phi) * np.sin(theta);
210+
# vy = vr * np.sin(phi) * np.cos(theta) + vphi * np.cos(phi)- vtheta * np.sin(phi) * np.sin(theta);
211+
# vz = vr * np.sin(theta) + vtheta * np.cos(theta);
212+
# return np.array([vx, vy, vz])
213+
214+
# # two more from jeff's code, plus my reverse of his
215+
# def jeff_transform_pos(ra, dec, dist,
216+
# R =np.array([[-0.05487395617553902, -0.8734371822248346,-0.48383503143198114],
217+
# [0.4941107627040048, -0.4448286178025452,0.7469819642829028],
218+
# [-0.8676654903323697, -0.1980782408317943,0.4559842183620723]]),
219+
# H =np.array([[0.9999967207734917, 0.0, 0.002560945579906427],
220+
# [0.0, 1.0, 0.0],
221+
# [-0.002560945579906427, 0.0, 0.9999967207734917]]),
222+
# offsett = np.array([8.112,0,0.0208])):
223+
# r_icrs = s2c_pos(dist,dec,ra)
224+
# r_gal = R @ r_icrs
225+
# r_gal -= offsett
226+
# r_gal = H @ r_gal
227+
# return r_gal
228+
# def jeff_reverse_transform_pos(r_gal,
229+
# R=np.array([[-0.05487395617553902, -0.8734371822248346,-0.48383503143198114],
230+
# [0.4941107627040048, -0.4448286178025452,0.7469819642829028],
231+
# [-0.8676654903323697, -0.1980782408317943,0.4559842183620723]]),
232+
# H=np.array([[0.9999967207734917, 0.0, 0.002560945579906427],
233+
# [0.0, 1.0, 0.0],
234+
# [-0.002560945579906427, 0.0, 0.9999967207734917]]),
235+
# offsett=np.array([8.112, 0, 0.0208]),return_cart=False):
236+
# r_helio = np.linalg.inv(H) @ r_gal
237+
# r_helio += offsett
238+
# r_icrs = np.linalg.inv(R) @ r_helio
239+
# if return_cart:
240+
# return r_icrs
241+
# else:
242+
# dist, dec, ra = c2s_pos(*r_icrs) # assumes `c2s_pos` converts Cartesian to spherical
243+
# return np.array([ra,dec,dist])
244+
245+
# def jeff_transform_vel(ra, dec, dist, pmra, pmdec, vlos,
246+
# R = np.array([[-0.05487395617553902, -0.8734371822248346,-0.48383503143198114],
247+
# [0.4941107627040048, -0.4448286178025452,0.7469819642829028],
248+
# [-0.8676654903323697, -0.1980782408317943,0.4559842183620723]]),
249+
# H = np.array([[0.9999967207734917, 0.0, 0.002560945579906427],
250+
# [0.0, 1.0, 0.0],
251+
# [-0.002560945579906427, 0.0, 0.9999967207734917]]),
252+
# offsett=np.array([8.112, 0, 0.0208]), solarmotion=np.array([12.9/100, 245.6/100, 7.78/100]),return_cart=True):
253+
# vlos /= 100
254+
# conversion_factor = 4.740470463533349
255+
# dist_with_conversion = dist * conversion_factor
256+
257+
# vra = dist_with_conversion * pmra / 100
258+
# vdec = dist_with_conversion * pmdec / 100
259+
260+
# r_gal = jeff_transform_pos(ra, dec, dist, R, H, offsett)
261+
# v_icrs = s2c_vel(dist, dec, ra, vlos, vdec, vra)
262+
# v_gal = H @ R @ v_icrs
263+
# v_gal += solarmotion
264+
# if return_cart:return v_gal
265+
# else:
266+
# v_sph = c2s_vel(r_gal[0], r_gal[1], r_gal[2],v_gal[0], v_gal[1], v_gal[2])
267+
# return v_sph
268+
# def jeff_reverse_transform_vel(r_gal, v_gal, R=np.array([[-0.05487395617553902, -0.8734371822248346,-0.48383503143198114],
269+
# [0.4941107627040048, -0.4448286178025452,0.7469819642829028],
270+
# [-0.8676654903323697, -0.1980782408317943,0.4559842183620723]]),
271+
# H=np.array([[0.9999967207734917, 0.0, 0.002560945579906427],
272+
# [0.0, 1.0, 0.0],
273+
# [-0.002560945579906427, 0.0, 0.9999967207734917]]),
274+
# offsett=np.array([8.112, 0, 0.0208]), solarmotion=np.array([12.9/100, 245.6/100, 7.78/100])):
275+
# v_gal -= solarmotion
276+
# v_icrs = np.linalg.inv(H @ R) @ v_gal
277+
# xyz_helio = jeff_reverse_transform_pos(r_gal,return_cart=True)
278+
# ra,dec,dist = jeff_reverse_transform_pos(r_gal,return_cart=False)
279+
# vlos, vdec, vra = c2s_vel(*xyz_helio, *v_icrs)
280+
# pmra = vra / (dist * 4.740470463533349) * 100
281+
# pmdec = vdec / (dist * 4.740470463533349) * 100
282+
# vlos *= 100
283+
# return np.array([vlos, pmra, pmdec])

‎df_sampling/core_imports.py‎

Lines changed: 41 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,41 @@
1+
# core_imports.py
2+
3+
# Core libraries
4+
import numpy as np
5+
import pandas as pd
6+
import matplotlib.pyplot as plt
7+
import seaborn as sns
8+
9+
sns.set_style('white')
10+
sns.set_palette('Dark2')
11+
12+
import astropy.units as u
13+
import astropy.constants as const
14+
from astropy import coordinates as Coords
15+
from astropy.coordinates import SkyCoord
16+
17+
from scipy.stats import linregress, pareto, uniform, halfcauchy, norm, multivariate_normal
18+
from scipy.special import gamma as gammaf
19+
from scipy.special import gammaln
20+
from scipy.optimize import fsolve, minimize
21+
from scipy.interpolate import interp1d
22+
from scipy.special import hyp2f1, betainc
23+
24+
from dataclasses import dataclass
25+
import glob as glob
26+
27+
import time
28+
import datetime
29+
30+
import sys
31+
import argparse
32+
# import cmdstanpy
33+
# import arviz as az
34+
import os
35+
import warnings
36+
import logging
37+
import re
38+
# import multiprocessing
39+
# import concurrent.futures
40+
# from concurrent.futures import ProcessPoolExecutor, as_completed
41+

0 commit comments

Comments
 (0)