Newer
Older
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
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
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
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
import os
import yaml
import numpy as np
import scipy.ndimage as nd
from astropy.io import fits
import matplotlib.pyplot as plt
import math
from CpicImgSim.config import cpism_refdata, solar_spectrum, MAG_SYSTEM
from CpicImgSim.utils import region_replace, random_seed_select
from CpicImgSim.io import log
from CpicImgSim.optics import filter_throughput
def sky_frame_maker(band, skybg, platescale, shape):
"""
generate a sky background frame.
Parameters
----------
band : str
The band of the sky background.
skybg : str
The sky background file name.
platescale : float
The platescale of the camera in arcsec/pixel.
shape : tuple
The shape of the output frame. (y, x)
Returns
-------
sky_bkg_frame : numpy.ndarray
The sky background frame.
"""
filter = filter_throughput(band)
sk_spec = solar_spectrum.renorm(skybg, MAG_SYSTEM, filter)
sky_bkg_frame = np.zeros(shape)
sky_bkg_frame += (sk_spec * filter).integrate() * platescale**2
return sky_bkg_frame
class CRobj(object):
"""
Cosmic ray object.
Attributes
----------
flux : float
The flux of the cosmic ray.
angle : float
The angle of the cosmic ray.
sigma : float
The width of the cosmic ray.
length : int
The length of the cosmic ray.
"""
def __init__(self, flux, angle, sigma, length) -> None:
self.flux = flux
self.angle = angle
self.sigma = sigma
self.length = length
class CosmicRayFrameMaker(object):
"""
Cosmic ray frame maker.
Parameters
----------
depth : float
The depth of the camera pixel in um.
pitch : float
The pitch of the camera pixel in um.
cr_rate : float
The cosmic ray rate per second per cm2.
"""
def __init__(self) -> None:
self.tmp_size = [7, 101]
self.freq_power = -0.9
self.trail_std = 0.1
self.depth = 10 # um
self.pitch = 13 # um
self.cr_rate = 1 # particle per s per cm2 from Miles et al. 2021
def make_CR(self, length, sigma, seed=-1):
"""
make a image of cosmic ray with given length and sigma.
Parameters
----------
length : int
The length of the cosmic ray in pixel.
sigma : float
The width of the cosmic ray in pixel.
Returns
-------
output : numpy.ndarray
The image of cosmic ray.
"""
h = self.tmp_size[0]
w = self.tmp_size[1]
freq = ((w-1)/2-np.abs(np.arange(w)-(w-1)/2)+1)**(self.freq_power)
x = np.arange(w) - (w-1)/2
hl = (length-1)/2
x_wing = np.exp(-(np.abs(x)-hl)**2/sigma/sigma/2)
x_wing[np.abs(x) < hl] = 1
cr = np.zeros([h, w])
center = (h-1)/2
for i in range(h):
phase = np.random.rand(w)*2*np.pi
trail0 = abs(np.fft.fft(freq*np.sin(phase) + 1j*x*np.cos(phase)))
# TODO maybe somthing wrong
trail_norm = (trail0 - trail0.mean())/trail0.std()
cr[i, :] = np.exp(-(i - center)**2/sigma/sigma/2) \
* (trail_norm * self.trail_std + 1) * x_wing
output = np.zeros([w, w])
d = (w-h)//2
output[d:d+h, :] = cr
return output
def _length_rand(self, N, seed=-1):
"""
randomly generate N cosmic ray length.
"""
len_out = []
seed = random_seed_select(seed=seed)
log.debug(f"cr length seed: {seed}")
for i in range(N):
x2y2 = 2
while x2y2 > 1:
lx = 1 - 2 * np.random.rand()
ly = 1 - 2 * np.random.rand()
x2y2 = lx * lx + ly * ly
z = 1 - 2 * x2y2
r = 2 * np.sqrt(x2y2 * (1 - x2y2))
length = abs(r / z * self.depth)
pitch = self.pitch
len_out.append(int(length/pitch))
return np.array(len_out)
def _number_rand(self, expt, pixsize, random=False, seed=-1):
"""
randomly generate the number of cosmic rays.
"""
area = (self.pitch / 1e4)**2 * pixsize[0] * pixsize[1]
ncr = area * expt * self.cr_rate
if random:
seed = random_seed_select(seed=seed)
log.debug(f"cr count: {seed}")
ncr = np.random.poisson(ncr)
else:
ncr = int(ncr)
self.ncr = ncr
return ncr
def _sigma_rand(self, N, seed=-1):
"""
randomly generate N cosmic ray sigma.
"""
sig_sig = 0.5 # asuming the sigma of size of cosmic ray is 0.5
seed = random_seed_select(seed=seed)
log.debug(f"cr width seed: {seed}")
sig = abs(np.random.randn(N))*sig_sig + 1/np.sqrt(12) * 1.2
# assume sigma is 1.2 times of pictch sig
return sig
def _flux_rand(self, N, seed=-1):
"""
randomly generate N cosmic ray flux.
"""
seed = random_seed_select(seed=seed)
log.debug(f"cr flux seed: {seed}")
u = np.random.rand(N)
S0 = 800
lam = 0.57
S = (-np.log(1-u)/lam + S0**0.25)**4
return S
def random_CR_parameter(self, expt, pixsize):
"""
randomly generate cosmic ray parameters, including number, length, flux, sigma and angle.
Parameters
----------
expt : float
The exposure time in second.
pixsize : list
The size of the image in pixel.
Loading
Loading full blame…