~ubuntu-branches/ubuntu/wily/eso-midas/wily-proposed

« back to all changes in this revision

Viewing changes to stdred/long/proc/lnhough.prg

  • Committer: Package Import Robot
  • Author(s): Ole Streicher
  • Date: 2014-04-22 14:44:58 UTC
  • Revision ID: package-import@ubuntu.com-20140422144458-okiwi1assxkkiz39
Tags: upstream-13.09pl1.2+dfsg
ImportĀ upstreamĀ versionĀ 13.09pl1.2+dfsg

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
! @(#)lnhough.prg       19.1 (ESO-DMD) 02/25/03 14:24:46
 
2
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
3
!.COPYRIGHT   (C) 1993 European Southern Observatory
 
4
!.IDENT       lnhough.prg
 
5
!.AUTHORS     Pascal Ballester (ESO/Garching)
 
6
!.KEYWORDS    Spectroscopy, Long-Slit
 
7
!.PURPOSE     Tutorial for context Long . PB 12.02.1993
 
8
!.VERSION     1.0  Command Creation  09-APR-1993
 
9
! 020906        last modif
 
10
!-------------------------------------------------------
 
11
!
 
12
CROSSREF wdisp wcent ystart line cat  mode  range
 
13
!
 
14
define/param p1  {avdisp}  ?    "Dispersion, tolerance, accuracy:"
 
15
define/param p2  {wcenter} ?    "Central wav., tolerance, accuracy:"
 
16
define/param p3  {ystart}  NUM    "Reference row:"
 
17
define/param p4  {lintab}  TAB    "Line table:"
 
18
define/param p5  {lincat}  TAB    "Line catalog:"
 
19
define/param p6   LINEAR   CHAR   "Mode:"
 
20
define/param p7   2.5      NUM    "Range:"
 
21
define/param p8   K        CHAR   "Keep/Visualise result (K/N/V)"
 
22
 
 
23
define/local  wcent/d/1/3  0.,0.,0.
 
24
define/local  wdisp/d/1/3  0.,0.,0.
 
25
define/local  werror/i/1/1  1
 
26
define/local  dim/i/1/1   0
 
27
 
 
28
write/keyw     wcent {P2}
 
29
write/keyw     wdisp {P1}
 
30
 
 
31
if wdisp(2) .eq. 0.  wdisp(2) = 15  ! Tolerance in percent
 
32
if wdisp(3) .eq. 0   wdisp(3) = 1   ! Accuracy in percent
 
33
if wcent(2) .eq. 0   wcent(2) = 50  ! Tolerance 100 pixels
 
34
if wcent(3) .eq. 0   wcent(3) = 0.2 ! Accuracy in pixels
 
35
 
 
36
set/format  f8.6
 
37
write/out "Initial values:"
 
38
write/out " Central Wav.  : {wcent(1)} +/- {wcent(2)} pix, {wcent(3)} pix"
 
39
write/out " Average disp. : {wdisp(1)} +/- {wdisp(2)} %, {wdisp(3)} %"
 
40
 
 
41
wdisp(2) = wdisp(2)/100 
 
42
wdisp(3) = wdisp(3)/100
 
43
!
 
44
define/local nblines/I/1/1 0
 
45
define/local conf/R/1/1 0.
 
46
!
 
47
@s lnident,seline
 
48
nblines = outputi(1)
 
49
!
 
50
copy/table {lintab} &l
 
51
select/table {lintab} ALL
 
52
create/column &l :WAVEC R*8 F10.3
 
53
 
 
54
 
 
55
define/local xc/d/1/1 0.
 
56
define/local xmean/r/1/1 0.
 
57
define/local range/r/1/1 0.
 
58
define/local factor/r/1/1 0.
 
59
define/local offset/r/1/1 0.
 
60
 
 
61
range  = {p7}*step(1)
 
62
 
 
63
! Dispersion relation is relative to central positions of the group
 
64
! of lines in pixel space.
 
65
statistics/table {lintab} :X
 
66
 
 
67
xmean = (outputr(3) - start(1))/step(1) + 1.
 
68
xc = npix(1)/2. + 0.5
 
69
 
 
70
factor = xmean
 
71
offset = (xmean - xc)*wdisp(1)
 
72
 
 
73
write/out "XCENTER = {xc}  XMEAN = {xmean} OFFSET={offset}"
 
74
 
 
75
wcent(1) = wcent(1) + offset
 
76
 
 
77
set/format f20.10
 
78
compute/table  &l  :XMEAN  =  (:X - ({start(1)}))/{step(1)}+1.-{factor}
 
79
compute/table  &l  :XNORM  =  :XMEAN/{factor}
 
80
 
 
81
define/local wlim/d/1/2  0.,0.
 
82
if lincat(1:1) .ne. "@" copy/table {lincat} &c
 
83
 
 
84
define/local hg_start/D/1/3 0.,0.,0.  ? +lower
 
85
define/local hg_step/D/1/3  0.,0.,0.  ? +lower
 
86
define/local hg_npix/I/1/3  0,0,0     ? +lower
 
87
 
 
88
if p6(1:1) .eq. "L"  then
 
89
  hg_start(1) = wdisp(1) * (1.-wdisp(2))
 
90
  hg_step(1)  = wdisp(1) * wdisp(3)
 
91
  hg_npix(1)  = 2.*wdisp(2)/wdisp(3)
 
92
 
 
93
  hg_start(2) = wcent(1) - wcent(2)*wdisp(1)
 
94
  hg_step(2)  = wcent(3)*wdisp(1)
 
95
  hg_npix(2)  = 2.*wcent(2)/wcent(3)
 
96
 
 
97
  dim = 2
 
98
  write/keyw inputc/c/1/10  :XMEAN
 
99
endif
 
100
 
 
101
if p6(1:1) .eq. "1" then
 
102
      dim = 1
 
103
      hg_start(1) = wcent(1) - wcent(2)*wdisp(1)
 
104
      hg_step(1)  = wcent(3)*wdisp(1)
 
105
      hg_npix(1)  = 2.*wcent(2)/wcent(3)
 
106
      inputr(2) = wdisp(1)
 
107
      write/keyw inputc/c/1/10 :XMEAN
 
108
endif
 
109
 
 
110
if p6(1:1) .eq. "N" then
 
111
      dim = 2
 
112
      inputr(2) = wdisp(1)*factor
 
113
      write/keyw inputc/c/1/10  :XNORM
 
114
 
 
115
      hg_start(1) = -0.30
 
116
      hg_step(1)  = m$abs(hg_start(1))/100.
 
117
      hg_npix(1)  = 200
 
118
 
 
119
      hg_start(2) = wcent(1) - wcent(2)*wdisp(1)
 
120
      hg_step(2)  = wcent(3)*wdisp(1)
 
121
      hg_npix(2)  = 2.*wcent(2)/wcent(3)
 
122
endif
 
123
 
 
124
if p6(1:1) .eq. "3" then
 
125
      dim = 3
 
126
      inputr(2) = wdisp(1)*factor
 
127
      range     = range/factor
 
128
      write/keyw inputc/c/1/10  :XNORM
 
129
 
 
130
      hg_start(1) = inputr(2) * (1.-wdisp(2))
 
131
      hg_step(1)  = inputr(2) * wdisp(3)
 
132
      hg_npix(1)  = 2.*wdisp(2)/wdisp(3)
 
133
 
 
134
      hg_start(2) = wcent(1) - wcent(2)*wdisp(1)
 
135
      hg_step(2)  = wcent(3)*wdisp(1)
 
136
      hg_npix(2)  = 2.*wcent(2)/wcent(3)
 
137
 
 
138
      hg_start(3) = -0.2
 
139
      hg_step(3)  = 0.01
 
140
      hg_npix(3)  = 40
 
141
endif
 
142
 
 
143
if dim .eq. 0  then
 
144
   write/out "Wrong Method : {p6}"
 
145
   return/exit
 
146
endif
 
147
 
 
148
write/keyw inputc/c/10/10  :WAVE
 
149
write/keyw in_a middumml.tbl
 
150
 
 
151
if lincat(1:1) .ne. "@" then
 
152
   write/keyw in_b   middummc.tbl
 
153
else
 
154
   write/keyw in_b   {LINCAT}
 
155
endif
 
156
if p8(1:1) .ne. "N" then
 
157
    write/keyw out_a  middummh.bdf
 
158
else
 
159
    write/keyw out_a @@@
 
160
endif
 
161
write/key out_b  {P6}
 
162
 
 
163
@s lnhough,bounds
 
164
 
 
165
inputi(4) = dim
 
166
inputr(4) = range
 
167
copy/keyw hg_npix/I/1/3  inputi/I/1/3
 
168
copy/keyw hg_start/D/1/3 inputd/D/1/3 
 
169
copy/keyw hg_step/D/1/3  INPUTD/D/4/3
 
170
 
 
171
run STD_EXE:lnhough.exe
 
172
 
 
173
define/local panmax/i/1/1 {outputi(5)}
 
174
 
 
175
conf = (outputr(1)/nblines)*100.
 
176
if conf .gt. 100. conf = 100.
 
177
 
 
178
if p8(1:1) .eq. "V"  then
 
179
   extract/image &y = &h[<,<,@{panmax}:>,>,@{panmax}]
 
180
   load/image &y scale=5,1
 
181
endif
 
182
 
 
183
define/local staep/D/1/3  0.,0.,0.
 
184
 
 
185
if p6(1:1) .eq. "1" then
 
186
   staep(1) = hg_start(1) + (outputi(3)-1)*hg_step(1)
 
187
   wcenter = staep(1)  - OFFSET
 
188
endif
 
189
 
 
190
if p6(1:1) .eq. "L" .or. p6(1:1) .eq. "N" then
 
191
   staep(1) = hg_start(1) + (outputi(3)-1)*hg_step(1)
 
192
   staep(2) = hg_start(2) + (outputi(4)-1)*hg_step(2)
 
193
 
 
194
   if p6(1:1) .eq. "L" then
 
195
       compute/table &l :WAVEC = {staep(2)}+{staep(1)}*:XMEAN
 
196
       avdisp   = staep(1)
 
197
       beta     = 0.
 
198
   endif
 
199
 
 
200
   if p6(1:1) .eq. "N" then
 
201
       compute/table &l :WAVEC = -
 
202
         {staep(2)}+{wdisp(1)}*{factor}*:XNORM+{staep(1)}*:XNORM*:XNORM
 
203
       beta     = staep(1)
 
204
   endif
 
205
 
 
206
   wcenter  = staep(2) - offset
 
207
endif
 
208
 
 
209
if p6(1:1) .eq. "3" then
 
210
   staep(1) = hg_start(1) + (outputi(3)-1)*hg_step(1)
 
211
   staep(2) = hg_start(2) + (outputi(4)-1)*hg_step(2)
 
212
   staep(3) = hg_start(3) + (outputi(5)-1)*hg_step(3)
 
213
 
 
214
   compute/table &L :WAVEC = -
 
215
     {staep(2)}+{staep(1)}*:XNORM+{staep(3)}*:XNORM*:XNORM
 
216
 
 
217
   avdisp   = staep(1)/factor
 
218
   wcenter  = staep(2) - offset
 
219
   beta     = staep(3)
 
220
endif
 
221
 
 
222
 
 
223
set/format f12.5
 
224
write/out "Confidence Level: {CONF} %"
 
225
write/out "Set values WCENTER = {WCENTER} and AVDISP = {AVDISP} and BETA={BETA}"
 
226
 
 
227
! Now Updates column :WAVEC in line.tbl
 
228
 
 
229
regres/poly  &l  :WAVEC  :X  2
 
230
copy/kk outputd/d/1/3  dispcoe/d/1/3
 
231
save/regress {lintab}  HREG
 
232
 
 
233
outputr(1) = conf
 
234
 
 
235
 
 
236
 
 
237
ENTRY BOUNDS
 
238
 
 
239
define/local dim/I/1/1 0
 
240
define/local hgs/D/1/1 0.
 
241
 
 
242
do dim = 1 3
 
243
   hgs = hg_step({dim})
 
244
   if hgs .lt. 0. then
 
245
      hg_start({dim}) = hg_start({dim}) + (hg_npix({dim})-1.)*hgs
 
246
      hg_step({dim}) = hgs*(-1.)
 
247
   endif
 
248
enddo
 
249