Mercurial > repos > astroteam > crbeam_astro_tool

annotate mc_rotate.py @ 0:f40d05521dca draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
author astroteam
date Fri, 25 Apr 2025 19:33:20 +0000
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
1 import math
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
2
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
3 import numpy as np
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
4
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
5
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
6 def mc_rotate(
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
7 mc_file, jet_half_size_deg, jet_direction_deg, psf_deg=180, rotated_file=None
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
8 ):
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
9 if rotated_file is None:
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
10 rotated_file = f"{mc_file}_{jet_half_size_deg}_{jet_direction_deg}_{psf_deg}"
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
11
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
12 jet_half_size = (
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
13 np.longdouble(jet_half_size_deg) / 180.0 * math.pi
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
14 ) # jet openning angle
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
15 jet_direction = (
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
16 np.longdouble(jet_direction_deg) / 180.0 * math.pi
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
17 ) # angle between direction to source and jet axes
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
18
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
19 if jet_half_size > math.pi or jet_half_size <= 0:
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
20 jet_half_size = math.pi
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
21 # here we assume that jet axis is located in y-z plane if angle=0, then jet points towards observer
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
22 jet_axes = np.array(
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
23 [0, np.sin(jet_direction), np.cos(jet_direction)], dtype=np.longdouble
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
24 )
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
25
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
26 if jet_direction > math.pi or jet_direction < 0:
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
27 raise ValueError(
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
28 "invalid jet_direction_deg param [0, 180] range expected")
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
29
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
30 with open(mc_file, "r") as lines, open(rotated_file, "w") as rotated:
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
31 column_names = []
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
32 for line in lines:
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
33 columns = line.split()
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
34 if len(columns) == 0:
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
35 continue
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
36 if columns[0] == "#":
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
37 if "weight" in columns:
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
38 column_names = (
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
39 columns[1:3]
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
40 + ["Theta", "Phi", "dT_raw/T", "dT_calc/T"]
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
41 + columns[9:]
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
42 )
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
43 print("# " + "\t".join(column_names), file=rotated)
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
44 else:
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
45 print(line.strip(), file=rotated)
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
46 else:
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
47 particle = np.array(list(map(np.longdouble, columns)))
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
48 beta = phi = 0 #dTcalc = 0
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
49
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
50 r = particle[2:5]
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
51 dz = r[
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
52 2
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
53 # dz is roughly equal to time delay in units of full path (assuming very small deflection angle)
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
54 ]
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
55
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
56 # negative dz may apper in data as a result of rounding errors
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
57 dz = max(dz, 0) # replacing negative values by zeros
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
58
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
59 r[2] = 1.0 - dz # now r is position of particle
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
60
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
61 # direction of momentum
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
62 n = particle[5:8]
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
63 n[2] = 1.0 - n[2]
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
64 # compensate rounding errors for n
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
65 mod_n = np.dot(n, n) ** 0.5
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
66 n = 1 / mod_n * n
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
67 t = 0 # calculated delay
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
68 if dz > 0: # there is some delay
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
69 r2 = np.dot(r, r)
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
70 # assert r2 < 1
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
71 if r2 > 1:
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
72 r *= 1.0 / np.sqrt(r2)
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
73 r2 = 1.0
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
74 # we continue trajectory in the direction of momentun until it reaches the sphere of radius 1.
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
75 # to find the delay t, we need to solve the following quadratic equation with respect to variable t:
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
76 # [(x + V_x)*t]^2 + [(y + V_y)*t]^2 + [(z + V_z)*t]^2 = 1
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
77 # with extra condition t > 0
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
78 # below is the solution:
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
79 rV = np.dot(n, r)
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
80 t = np.sqrt(1.0 - r2 + rV * rV) - rV
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
81
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
82 # calculating the end point of the trajectory
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
83 r = r + t * n
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
84
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
85 # normalizing r to 1 to fix the rounding errors
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
86 modR = np.dot(r, r) ** 0.5
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
87 r = 1.0 / modR * r
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
88
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
89 # now we will construct the rotation of coordinates which converts vector r to 0,0,1
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
90 # to find the starting direction of the particle momentum
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
91 dr2 = r[0] ** 2 + r[1] ** 2
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
92 if dr2 > 0:
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
93 # building rotation matrix to move x,y,z to 0,0,1
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
94 norm = 1.0 / (dr2**0.5)
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
95 # rotation axes:
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
96 ux = norm * r[1]
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
97 uy = -norm * r[0]
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
98 # uz=0
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
99 cos = r[2]
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
100 sin = (1 - cos**2) ** 0.5
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
101
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
102 toCenter = np.array(
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
103 [
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
104 [cos + ux * ux * (1 - cos), ux *
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
105 uy * (1 - cos), uy * sin],
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
106 [ux * uy * (1 - cos), cos + uy *
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
107 uy * (1 - cos), -ux * sin],
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
108 [-uy * sin, ux * sin, cos],
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
109 ]
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
110 )
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
111
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
112 zAxes = np.array([0.0, 0.0, 1.0], dtype=np.longdouble)
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
113 # test
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
114 error = np.fabs(np.dot(toCenter, r) - zAxes).max()
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
115 assert error < 1e-7
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
116 # finding starting direction of the momentum
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
117 nInitial = np.dot(toCenter, zAxes)
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
118 norm_nInitial = np.sqrt(np.dot(nInitial, nInitial))
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
119 nInitial *= (
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
120 1.0 / norm_nInitial
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
121 ) # normalize after rotation to avoid rounding errors
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
122 # calculating angle between starting direction and jet axis
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
123 startCos = np.dot(nInitial, jet_axes)
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
124 startAngle = 0
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
125 if (
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
126 startCos < 1
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
127 ): # avoid roundout errors if angle is very close to zero
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
128 startAngle = np.arccos(startCos)
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
129 # TODO: here instead of using step function we can modify weight depending
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
130 # on angular distance from jet axes
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
131 if (
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
132 startAngle < jet_half_size
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
133 ): # if the particle is inside the cone (passes cut)
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
134 # calculate angle between arrival direction and direction to source
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
135 cosBeta = np.dot(
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
136 r, n
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
137 ) # this angle is invariant (doesn't depend on turns)
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
138 if (
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
139 cosBeta < 1
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
140 ): # avoid rounding errors if angle is very close to zero
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
141 beta = np.arccos(cosBeta)
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
142 nFinal = np.dot(toCenter, n).transpose()
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
143 cosPhi = nFinal[0] / np.sin(
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
144 beta
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
145 # second angle (could be useful in assymetric tasks) if beta=0 it is undefined (could be any)
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
146 )
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
147 if cosPhi < 1 and cosPhi > -1:
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
148 phi = np.arccos(
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
149 cosPhi
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
150 ) # returns value in interval [0,pi]
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
151 if nFinal[1] < 0:
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
152 phi = 2 * math.pi - phi
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
153
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
154 # s = math.sin(startAngle+beta)
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
155 # if s>0.:
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
156 # dTcalc = (math.sin(startAngle)+math.sin(beta))/s - 1.
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
157 # dTcalc = (1.-min(modR, 1.))/cosBeta # condition modR > 1 can only occur due to rounding error
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
158 else:
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
159 continue
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
160 elif (
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
161 jet_direction > jet_half_size
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
162 ): # no deflection: check if jet is pointing towards observer
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
163 continue
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
164
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
165 if beta * 180.0 / math.pi > psf_deg:
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
166 continue
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
167
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
168 output = [
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
169 particle[0],
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
170 particle[1],
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
171 beta * 180.0 / math.pi,
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
172 phi * 180.0 / math.pi,
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
173 dz,
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
174 t,
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
175 ] + list(particle[8:])
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
176 output = ["{0:g}".format(i) for i in output]
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
177 print(*output, file=rotated)
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
178
f40d05521dca planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
astroteam
parents:
diff changeset
179 return rotated_file