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
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
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
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
|
C------------------------------------------------------------------------------
subroutine ftwldp(xpix,ypix,xref,yref,xrefpix,yrefpix,
& xinc,yinc,rot,type,xpos,ypos,status)
C Fortran version of worldpos.c -- WCS Algorithms from Classic AIPS
C Translated by James Kent Blackburn -- HEASARC/GSFC/NASA -- November 1994
C routine to determine accurate position from pixel coordinates
C does: -SIN, -TAN, -ARC, -NCP, -GLS, -MER, -AIT projections
C returns 0 = good,
C 501 = angle too large for projection;
C Input:
C dbl xpix x pixel number (RA or long without rotation)
C dbl ypiy y pixel number (dec or lat without rotation)
C dbl xref x reference coordinate value (deg)
C dbl yref y reference coordinate value (deg)
C dbl xrefpix x reference pixel
C dbl yrefpix y reference pixel
C dbl xinc x coordinate increment (deg)
C dbl yinc y coordinate increment (deg)
C dbl rot rotation (deg) (from N through E)
C chr type projection type code e.g. "-SIN"
C Output:
C dbl xpos x (RA) coordinate (deg)
C dbl ypos y (dec) coordinate (deg)
C int status error status flag, zero
integer status
double precision xpix,ypix,xref,yref,xrefpix,yrefpix
double precision xinc,yinc,rot,xpos,ypos
character*(*) type
integer error1,error4
parameter (error1=501)
parameter (error4=504)
double precision cosr,sinr,dx,dy,dz,temp
double precision sins,coss,dect,rat,dt,l,m,mg,da,dd,cos0,sin0
double precision dec0,ra0,decout,raout
double precision geo1,geo2,geo3
double precision cond2r
parameter (cond2r=1.745329252d-2)
double precision twopi,deps
parameter (twopi = 6.28318530717959)
parameter (deps = 1.0d-5)
integer i,itype
character*4 ctypes(8)
data ctypes/ '-SIN', '-TAN', '-ARC', '-NCP',
& '-GLS', '-MER', '-AIT', '-STG' /
if (status .gt. 0) return
C *** Offset from ref pixel
dx = (xpix-xrefpix) * xinc
dy = (ypix-yrefpix) * yinc
C *** Take out rotation
cosr = dcos(rot*cond2r)
sinr = dsin(rot*cond2r)
if (rot .ne. 0.0) then
temp = dx * cosr - dy * sinr
dy = dy * cosr + dx * sinr
dx = temp
end if
C *** Find type of coordinate transformation (0 is linear)
itype = 0
do 10 i = 1, 8
if (ctypes(i) .eq. type) itype = i
10 continue
C *** default, linear result for error return
xpos = xref + dx
ypos = yref + dy
C *** Convert to radians
ra0 = xref * cond2r
dec0 = yref * cond2r
l = dx * cond2r
m = dy * cond2r
sins = l*l + m*m
decout = 0.0
raout = 0.0
cos0 = dcos(dec0)
sin0 = dsin(dec0)
C *** Process by case
if (itype .eq. 0) then
C *** LINEAR
rat = ra0 + l
dect = dec0 + m
else if (itype .eq. 1) then
C *** SINE from '-SIN' type
if (sins .gt. 1.0) then
status = error1
goto 30
end if
coss = dsqrt(1.0 - sins)
dt = sin0 * coss + cos0 * m
if ((dt .gt. 1.0) .or. (dt .lt. -1.0)) then
status = error1
goto 30
end if
dect = dasin(dt)
rat = cos0 * coss - sin0 * m
if ((rat .eq. 0.0) .and. (l .eq. 0.0)) then
status = error1
goto 30
end if
rat = datan2 (l, rat) + ra0
else if (itype .eq. 2) then
C *** TANGENT from '-TAN' type
if (sins .gt. 1.0) then
status = error1
goto 30
end if
dect = cos0 - m * sin0
if (dect .eq. 0.0) then
status = error1
goto 30
end if
rat = ra0 + datan2(l, dect)
dect = datan(dcos(rat-ra0) * (m * cos0 + sin0) / dect)
else if (itype .eq. 3) then
C *** Arc from '-ARC' type
if (sins .ge. twopi * twopi / 4.0) then
status = error1
goto 30
end if
sins = dsqrt(sins)
coss = dcos(sins)
if (sins .ne. 0.0) then
sins = dsin(sins) / sins
else
sins = 1.0
end if
dt = m * cos0 * sins + sin0 * coss
if ((dt .gt. 1.0) .or. (dt .lt. -1.0)) then
status = error1
goto 30
end if
dect = dasin(dt)
da = coss - dt * sin0
dt = l * sins * cos0
if ((da .eq. 0.0) .and. (dt .eq. 0.0)) then
status = error1
goto 30
end if
rat = ra0 + datan2(dt, da)
else if (itype .eq. 4) then
C *** North Celestial Pole from '-NCP' type
dect = cos0 - m * sin0
if (dect .eq. 0.0) then
status = error1
goto 30
end if
rat = ra0 + datan2(l, dect)
dt = dcos(rat-ra0)
if (dt .eq. 0.0) then
status = error1
goto 30
end if
dect = dect / dt
if ((dect .gt. 1.0) .or. (dect .lt. -1.0)) then
status = error1
goto 30
end if
dect = dacos(dect)
if (dec0 .lt. 0.0) dect = -dect
else if (itype .eq. 5) then
C *** Global Sinusoid from '-GLS' type
dect = dec0 + m
if (dabs(dect) .gt. twopi/4.0) then
status = error1
goto 30
end if
coss = dcos(dect)
if (dabs(l) .gt. twopi*coss/2.0) then
status = error1
goto 30
end if
rat = ra0
if (coss .gt. deps) rat = rat + l / coss
else if (itype .eq. 6) then
C *** Mercator from '-MER' type
dt = yinc * cosr + xinc * sinr
if (dt .eq. 0.0) dt = 1.0
dy = (yref/2.0 + 45.0) * cond2r
dx = dy + dt / 2.0 * cond2r
dy = dlog(dtan(dy))
dx = dlog(dtan(dx))
geo2 = dt * cond2r / (dx - dy)
geo3 = geo2 * dy
geo1 = dcos(yref * cond2r)
if (geo1 .le. 0.0) geo1 = 1.0
rat = l / geo1 + ra0
if (dabs(rat - ra0) .gt. twopi) then
status = error1
goto 30
end if
dt = 0.0
if (geo2 .ne. 0.0) dt = (m + geo3) / geo2
dt = dexp(dt)
dect = 2.0 * datan(dt) - twopi / 4.0
else if (itype .eq. 7) then
C *** Aitoff from '-AIT' type
dt = yinc * cosr + xinc * sinr
if (dt .eq. 0.0) dt = 1.0
dt = dt * cond2r
dy = yref * cond2r
dx = dsin(dy+dt)/dsqrt((1.0+dcos(dy+dt))/2.0) -
& dsin(dy)/dsqrt((1.0+dcos(dy))/2.0)
if (dx .eq. 0.0) dx = 1.0
geo2 = dt / dx
dt = xinc * cosr - yinc * sinr
if (dt .eq. 0.0) dt = 1.0
dt = dt * cond2r
dx = 2.0 * dcos(dy) * dsin(dt/2.0)
if (dx .eq. 0.0) dx = 1.0
geo1 = dt * dsqrt((1.0+dcos(dy)*dcos(dt/2.0))/2.0) / dx
geo3 = geo2 * dsin(dy) / dsqrt((1.0+dcos(dy))/2.0)
rat = ra0
dect = dec0
if ((l .eq. 0.0) .and. (m .eq. 0.0)) goto 20
dz = 4.0-l*l/(4.0*geo1*geo1)-((m+geo3)/geo2)*((m+geo3)/geo2)
if ((dz .gt. 4.0) .or. (dz .lt. 2.0)) then
status = error1
goto 30
end if
dz = 0.5 * dsqrt(dz)
dd = (m+geo3) * dz / geo2
if (dabs(dd) .gt. 1.0) then
status = error1
goto 30
end if
dd = dasin(dd)
if (dabs(dcos(dd)) .lt. deps) then
status = error1
goto 30
end if
da = l * dz / (2.0 * geo1 * dcos(dd))
if (dabs(da) .gt. 1.0) then
status = error1
goto 30
end if
da = dasin(da)
rat = ra0 + 2.0 * da
dect = dd
else if (itype .eq. 8) then
C *** Stereographic from '-STG' type
dz = (4.0 - sins) / (4.0 + sins)
if (dabs(dz) .gt. 1.0) then
status = error1
goto 30
end if
dect = dz * sin0 + m * cos0 * (1.0+dz) / 2.0
if (dabs(dect) .gt. 1.0) then
status = error1
goto 30
end if
dect = dasin(dect)
rat = dcos(dect)
if (dabs(rat) .lt. deps) then
status = error1
goto 30
end if
rat = l * (1.0+dz) / (2.0 * rat)
if (dabs(rat) .gt. 1.0) then
status = error1
goto 30
end if
rat = dasin(rat)
mg = 1.0 + dsin(dect)*sin0 + dcos(dect)*cos0*dcos(rat)
if (dabs(mg) .lt. deps) then
status = error1
goto 30
end if
mg = 2.0 * (dsin(dect)*cos0 - dcos(dect)*sin0*dcos(rat)) / mg
if (dabs(mg-m) .gt. deps) rat = twopi/2.0 - rat
rat = ra0 + rat
else
C *** Unsupported Projection
status = error4
goto 30
end if
20 continue
C *** Return RA in range
raout = rat
decout = dect
if (raout-ra0 .gt. twopi/2.0) raout = raout - twopi
if (raout-ra0 .lt. -twopi/2.0) raout = raout + twopi
if (raout .lt. 0.0) raout = raout + twopi
C *** Correct units back to degrees
xpos = raout / cond2r
ypos = decout / cond2r
30 continue
end
|
