summaryrefslogtreecommitdiffstats
path: root/ast/ast_tester/testspecframe.f
blob: 4af260695695ea6e7a634ec28eb36c0aa4a62b3e (plain)
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
      program testspecframe
      implicit none
      include 'SAE_PAR'
      include 'AST_PAR'

      double precision rf, x, y
      integer status, sf, sf1, sf2, fs
      status = sai__ok

      sf = ast_specframe( 'system=freq,unit=Hz', status )
      if( ast_GetD( sf, 'SpecOrigin', status ) .ne. 0.0 ) then
         call stopit( status, 'Error 0' )
      end if

      rf = ast_GetD( sf, 'RestFreq', status )
      call ast_SetD( sf, 'SpecOrigin', rf*1.0D9, status )
      if( abs( ast_GetD( sf, 'SpecOrigin', status ) - rf*1.0D9 )
     :    .gt. 0.1 ) then
         call stopit( status, 'Error 1' )
      end if

      call ast_setc( sf, 'Unit(1)', 'GHz', status )
      if( ast_GetD( sf, 'SpecOrigin', status ) .ne. rf ) then
         call stopit( status, 'Error 2' )
      end if

      call checkdump( sf, 'Error 3', status )

      call ast_setc( sf, 'System', 'vrad', status )
      if( abs( ast_GetD( sf, 'SpecOrigin', status ) ) .gt. 1.0D-8 ) then
         write(*,*) ast_GetD( sf, 'SpecOrigin', status )
         call stopit( status, 'Error 4' )
      end if

      call ast_setc( sf, 'System', 'freq', status )
      call ast_setc( sf, 'Unit(1)', 'Hz', status )

      if( abs( ast_GetD( sf, 'SpecOrigin', status ) - rf*1.0D9 )
     :    .gt. 0.1 ) then
         write(*,*) ast_GetD( sf, 'SpecOrigin', status )
         call stopit( status, 'Error 5' )
      end if

      call ast_setc( sf, 'StdOfRest', 'LSRD', status )
      if( abs( ast_GetD( sf, 'SpecOrigin', status ) -
     :                  rf*1.00000212890848D9 ) .gt. 10.0 ) then
         write(*,*) ast_GetD( sf, 'SpecOrigin', status )
         write(*,*) 'Should be ',rf*1.00000212890848D9
         call stopit( status, 'Error 6' )
      end if


      sf1 = ast_specframe( 'system=freq,unit=Hz', status )
      call ast_setd( sf1, 'SpecOrigin', 1.0D20, status )
      sf2 = ast_specframe( 'system=freq,unit=Hz', status )
      call ast_setd( sf2, 'SpecOrigin', 1.01D20, status )
      fs = ast_convert( sf1, sf2, "", status );

      x = 0.03D20
      call ast_tran1( fs, 1, x, .true., y, status )
      if( abs( y - 0.02D20 ) .gt. 0.0 ) then
         write(*,*) y, y - 0.02D20
         call stopit( status, 'Error 7' )
      end if

      if( ast_getl( sf1, 'AlignSpecOffset', status ) ) then
         call stopit( status, 'Error 8' )
      end if
      call ast_setl( sf1, 'AlignSpecOffset', .true., status )
      call ast_setl( sf2, 'AlignSpecOffset', .true., status )

      fs = ast_convert( sf1, sf2, "", status );

      x = 0.03D20
      call ast_tran1( fs, 1, x, .true., y, status )
      if( abs( y - x ) .gt. 0.0 ) then
         write(*,*) y, y - x
         call stopit( status, 'Error 9' )
      end if

      sf = ast_specframe( 'system=freq,unit=Hz', status )
      call ast_setc( sf, 'SourceVRF', 'LSRK', status )
      call ast_setd( sf, 'SourceVel', 1000.0D0, status )

      call ast_setc( sf, 'SourceVRF', 'BARY', status )
      call ast_setc( sf, 'SourceSys', 'ZOPT', status )

      if( abs( ast_getd(  sf, 'SourceVel', status ) -
     :         0.00334028336870307D0 ) .gt. 1.0D-10 ) then
         write(*,*) ast_getd(  sf, 'SourceVel', status )
         call stopit( status, 'Error 11' )
      end if

      call checkdump( sf, 'Error 10', status )
      call ast_setc( sf, 'SourceVRF', 'LSRK', status )
      call ast_setc( sf, 'SourceSys', 'VREL', status )

      if( abs( ast_getd(  sf, 'SourceVel', status ) -
     :         1000.0D0 ) .gt. 1.0D-6 ) then
         write(*,*) ast_getd(  sf, 'SourceVel', status )
         call stopit( status, 'Error 12' )
      end if

      if( status .eq. sai__ok ) then
         write(*,*) 'All SpecFrame tests passed'
      else
         write(*,*) 'SpecFrame tests failed'
      end if

      end

      subroutine stopit( status, text )
      implicit none
      include 'SAE_PAR'
      integer status
      character text*(*)

      if( status .ne. sai__ok ) return
      status = sai__error
      write(*,*) text

      end


      subroutine checkdump( obj, text, status )
      implicit none
      include 'SAE_PAR'
      include 'AST_PAR'
      character text*(*)
      integer obj, status, next, end, ch, result, ll, overlap
      external mysource, mysink
      character buf*25000

      common /ss1/ buf
      common /ss2/ next, end, ll

      if( status .ne. sai__ok ) return

      ch = ast_channel( mysource, mysink, ' ', status )


      ll = 110
      next = 1
      if( ast_write( ch, obj, status ) .ne.1 ) then
         write(*,*) text
         call stopit( status, 'Cannot write supplied object to '//
     :                'channel' )
      end if

      next = 1
      result = ast_read( ch, status )
      if( result .eq. ast__null ) then
         write(*,*) text
         call stopit( status, 'Cannot read object from channel' )
      end if



      if( ast_getd( obj, 'specorigin', status ) .ne.
     :    ast_getd( result, 'specorigin', status ) ) then
         call ast_Show( obj, status )
         call ast_Show( result, status )
         write(*,*) text
         call stopit( status, 'Object has changed' )
      end if

      end

      subroutine sink1( status )
      implicit none
      include 'SAE_PAR'
      include 'AST_PAR'

      logical fsfound, done
      common /sink1com/ fsfound, done

      integer status, l
      character line*200

      if( status .ne. sai__ok ) return
      call ast_getline( line, l, status )

      if( index( line( : l ),'Unc =' ) .GT. 0 ) then
         done = .true.

      else if( .not. done .and.
     :         index( line( : l ),'FrameSet' ) .GT. 0 ) then
         fsfound= .true.
      end if

      end

      subroutine mysource( status )
      implicit none
      include 'SAE_PAR'
      include 'AST_PAR'
      integer status, next, end, ll
      character buf*25000

      common /ss1/ buf
      common /ss2/ next, end, ll

      if( status .ne. sai__ok ) return

      if( next .ge. end ) then
         call ast_putline( buf, -1, status )
      else
         call ast_putline( buf( next : ), ll, status )
      endif

      next = next + ll

      end

      subroutine mysink( status )
      implicit none
      include 'SAE_PAR'
      include 'AST_PAR'
      integer status, next, end, f, l, ll
      character buf*25000
      character line*1000

      common /ss1/ buf
      common /ss2/ next, end, ll

      if( status .ne. sai__ok ) return

      line = ' '
      call ast_getline( line, l, status )
      call chr_fandl( line( : l ), f, l )
      buf( next : ) = line( f : l )
      l = l - f + 1

      if( next + ll - 1 .ge. 25000 ) then
         write(*,*)
         call stopit( status, 'Buffer overflow in mysink!!' )
      else if( l .gt. ll ) then
         write(*,*)
         write(*,*) buf( next : next + l)
         write(*,*) 'Line length ',l
         call stopit( status, 'Line overflow in mysink!!' )
      else
         end = next + l
         buf( end : next + ll - 1 ) = ' '
      endif

      next = next + ll

      end