-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathsn1987a_model_lobes_ptsrc.pro
executable file
·147 lines (121 loc) · 4.22 KB
/
sn1987a_model_lobes_ptsrc.pro
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
function sn1987a_model_lobes_ptsrc,x,y,a,doplot=doplot
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;+
; This new model needs to be simplified and clearly define the characteristic radius.
; It should consist of a 2-D Gaussian elliptical (or circular?) torus
; plus 4 tangentially extended lobes with varying weights + central pointsource.
;
;
; Parameters:
; xpos = 1D array or x values
; ypos = 1D array of y values
; x0 = X coordinate of center of remnant
; y0 = Y coordinate of center of remnant
; sigr = Gaussian width of circular torus
; r0 = radius of peak brightness of torus
; a0 = background level of torus component
; a1 = weighting of torus component
; s0 = intensity of spot 1
; theta0 = position angle of spot 1
; sigt0 = tangential width of spot 1
; sigr0 = radial width of spot 1
; s1 = intensity of spot 2
; theta1 = position angle of spot 2
; sigt1 = tangential width of spot 2
; sigr1 = radial width of spot 2
; s2 = intensity of spot 3
; theta2 = position angle of spot 3
; sigt2 = tangential width of spot 3
; sigr2 = radial width of spot 3
; s3 = intensity of spot 4
; theta3 = position angle of spot 4
; sigt3 = tangential width of spot 4
; sigr3 = radial width of spot 4
; xsrc = location pointsource (tie to x0)
; ysrc = location pointsource (tie to y0)
; wsrc = width of the pointsource (in pixels)
; fsrc = flux of the source
;-
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
a0=a[0] ;background component
a1=a[1] ;scaling of torus
x0=a[2] ;x center
y0=a[3] ;y center
sigr=a[4] ;gaussian width of torus
r0=a[5] ;radius of peak brightness
s0=a[6] ;intensity of spot 1
theta0=a[7] ;position angle of spot 1
sigt0=a[8] ;tangential width of spot 1
sigr0=a[9] ;radial width of spot 1
s1=a[10] ;intensity of spot 2
theta1=a[11] ;position angle of spot 2
sigt1=a[12] ;tangential width of spot 2
sigr1=a[13] ;radial width of spot 2
s2=a[14] ;intensity of spot 3
theta2=a[15] ;position angle of spot 3
sigt2=a[16] ;tangential width of spot 3
sigr2=a[17] ;radial width of spot 3
s3=a[18] ;intensity of spot 4
theta3=a[19] ;position angle of spot 4
sigt3=a[20] ;tangential width of spot 4
sigr3=a[21] ;radial width of spot 4
counts=a[22]
xsrc=a[23]
ysrc=a[24]
wsrc=a[25]
fsrc=a[26]
; print, 'HERE: r0 = ',r0
;;Step 1: build Gaussian ellipse
xpos = (float(indgen(200)) - 100.0) - 0.5
xpos = rebin(xpos,200,200,/sample)
ypos = transpose(xpos)
r2=(xpos-x0)^2+(ypos-y0)^2
r=sqrt(r2)
; g=1./(2.*!pi*sigr^2)*exp(-(r-r0)^2/(2.*sigr^2))
g=exp(-(r-r0)^2/(2.*sigr^2))
; g=exp(-0.5*((r-r0)/sigr)^2)
image1=a0+a1*g
;;Step 2: add lobes
;each lobe on torus has position angle,intensity,tangential width (& radial width?)
;theta0,s0,sigt0,sigr0
theta = atan(ypos,xpos)
w=where(theta lt 0)
theta[w]=theta[w]+2*!pi
t0=theta-theta0
w=where(t0 gt !pi,nw)
if nw gt 0 then t0[w]=t0[w]-2.*!pi
t1=theta-theta1
w=where(t1 gt !pi,nw)
if nw gt 0 then t1[w]=t1[w]-2.*!pi
t2=theta-theta2
; w=where(t2 lt -!pi)
; t2[w]=t2[w]+2.*!pi
t3=theta-theta3
w=where(t3 lt -!pi,nw)
if nw gt 0 then t3[w]=t3[w]+2.*!pi
rsrc = sqrt((xpos - xsrc)^2d + (ypos - ysrc)^2d)
imsrc = fsrc[0]*exp(-rsrc^2d/(2d*wsrc^2d))
image2=s0*exp(-(t0)^2/(2.*sigt0^2)) ;*exp(-(r-r0)^2/(2.*sigr0^2imsrc = fsrc[0]*exp(-rsrc^2d/(2d*wsrc^2d)))) ;spot 1
image3=s1*exp(-(t1)^2/(2.*sigt1^2));*exp(-(r-r0)^2/(2.*sigr1^2)) ;spot 2
image4=s2*exp(-(t2)^2/(2.*sigt2^2));*exp(-(r-r0)^2/(2.*sigr2^2)) ;spot 3
image5=s3*exp(-(t3)^2/(2.*sigt3^2));*exp(-(r-r0)^2/(2.*sigr3^2)) ;spot 4
image6=image1*[image2+image3+image4+image5]+image1+imsrc ;image1 = ring, multiply to set lobes = 0 in center.
;;Step 4: Scale by total counts in image
image7=image6*counts/total(image6)
; doplot=1
if keyword_set(doplot) then begin
!p.multi=[0,3,3]
rdis,image1,/silent
rdis,image2,/silent
rdis,image3,/silent
rdis,image4,/silent
rdis,image5,/silent
; rdis,image6,/silent
rdis,imsrc,/silent
rdis,image7,/silent
!p.multi=0
endif
image=image7
;; stop
return,image
end