← Back to branch summary

~ubuntu-branches/debian/jessie/eso-midas/jessie

« back to all changes in this revision

Viewing changes to contrib/pepsys/src/meanstar.for

  • 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:
contrib/pepsys/src/meanstar.for
 
1
C===========================================================================
 
2
C Copyright (C) 1995-2010 European Southern Observatory (ESO)
 
3
C
 
4
C This program is free software; you can redistribute it and/or 
 
5
C modify it under the terms of the GNU General Public License as 
 
6
C published by the Free Software Foundation; either version 2 of 
 
7
C the License, or (at your option) any later version.
 
8
C
 
9
C This program is distributed in the hope that it will be useful,
 
10
C but WITHOUT ANY WARRANTY; without even the implied warranty of
 
11
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
12
C GNU General Public License for more details.
 
13
C
 
14
C You should have received a copy of the GNU General Public 
 
15
C License along with this program; if not, write to the Free 
 
16
C Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, 
 
17
C MA 02139, USA.
 
18
C
 
19
C Correspondence concerning ESO-MIDAS should be addressed as follows:
 
20
C       Internet e-mail: midas@eso.org
 
21
C       Postal address: European Southern Observatory
 
22
C                       Data Management Division 
 
23
C                       Karl-Schwarzschild-Strasse 2
 
24
C                       D 85748 Garching bei Muenchen 
 
25
C                       GERMANY
 
26
C===========================================================================
 
27
C
 
28
         PROGRAM MEANSTAR
 
29
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
30
C.IDENT         meanstar.for
 
31
C.MODULE
 
32
C.AUTHOR        Andrew T. Young
 
33
C.KEYWORD
 
34
C.LANGUAGE      FORTRAN 77
 
35
C.PURPOSE       collapse star table with multiple positions to means.
 
36
C.COMMENTS
 
37
C.RETURNS
 
38
C.ENVIRONMENT
 
39
C.VERSION       4.3
 
40
C
 
41
C 100618        last modif
 
42
C-----------------------------------------------------------------------------
 
43
C*****************************************************************************
 
44
C
 
45
C
 
46
C       Reads a sorted star table with duplicate positions, and tries to
 
47
C       estimate a "best" position for each star.  Complains if the spread
 
48
C       in either coordinate exceeds a minute of arc.
 
49
C
 
50
C       This is done in the command CONVERT/PHOT.
 
51
C
 
52
C       This program is modified from the "esodcon" program, and may contain
 
53
C       fossils from it.
 
54
C
 
55
C
 
56
C*****************************************************************************
 
57
C
 
58
C
 
59
        IMPLICIT NONE
 
60
C
 
61
C  BEGIN Declarations:
 
62
C
 
63
C
 
64
        INTEGER ITBL, NEWTBL
 
65
        INTEGER NCOLS,NROWS,NSORTC,NWPRAL,NROWSAL, ISTAT
 
66
C
 
67
C
 
68
        CHARACTER CARD*78, TBLFIL*80
 
69
        CHARACTER*32 OBJECT,OLDOBJ,CTEST
 
70
C
 
71
        INTEGER ITEM,ITEMS
 
72
        INTEGER KOBJ,KRA,KDEC,KEQ,NROW
 
73
        INTEGER JOBJ,JRA,JDEC,JEQ,NROUT
 
74
C
 
75
        INTEGER MXITEM
 
76
        PARAMETER (MXITEM=900)
 
77
        REAL RA,DEC,EQUINOX
 
78
        REAL RAS(MXITEM),DECS(MXITEM),EQUINOXS(MXITEM)
 
79
C
 
80
        LOGICAL NULL
 
81
C
 
82
C    Types for external fcns.:
 
83
C
 
84
C
 
85
C  Set up MIDAS declarations:
 
86
C
 
87
        INTEGER MADRID(1)
 
88
C
 
89
        COMMON /VMR/ MADRID
 
90
C
 
91
        INCLUDE 'MID_INCLUDE:ST_DEF.INC'
 
92
C
 
93
        INTEGER IUNIT
 
94
        INTEGER NACTEL,NULLS
 
95
C
 
96
C  END Declarations.
 
97
C
 
98
C
 
99
C  BEGIN DATA statements:
 
100
C
 
101
        INCLUDE 'MID_INCLUDE:ST_DAT.INC'
 
102
C
 
103
C
 
104
C  END DATA statements.
 
105
C
 
106
C
 
107
C  ********************  PROLOGUE  ********************
 
108
C
 
109
      CALL STSPRO ('MEANSTAR')
 
110
C
 
111
      OLDOBJ=' '
 
112
      NROW=0
 
113
      NROUT=0
 
114
C
 
115
C       Set up INPUT table file:
 
116
C
 
117
      CALL TV('Opening  sdata.tbl')
 
118
      CALL TBTOPN('sdata.tbl', 1,  ITBL, ISTAT)
 
119
      IF (ISTAT.NE.0) CALL TERROR (ITBL,1,'Could not open "sdata.tbl".')
 
120
C
 
121
      CALL TBIGET(ITBL,  NCOLS,NROWS,NSORTC,NWPRAL,NROWSAL, ISTAT)
 
122
      IF (ISTAT.NE.0) CALL TERROR
 
123
     1      (ITBL,1,'Could not get basic table data.')
 
124
C
 
125
C       Get column pointers...
 
126
C
 
127
      CALL TBLSER (ITBL, 'OBJECT',  KOBJ,ISTAT)
 
128
      IF (ISTAT.NE.0 .OR. KOBJ.EQ.-1)
 
129
     1    CALL TERROR(ITBL,2,'Could not find column OBJECT')
 
130
      CALL TBLSER (ITBL, 'RA',  KRA,ISTAT)
 
131
      IF (ISTAT.NE.0 .OR. KRA.EQ.-1)
 
132
     1    CALL TERROR(ITBL,3,'Could not find column RA')
 
133
      CALL TBLSER (ITBL, 'DEC',  KDEC,ISTAT)
 
134
      IF (ISTAT.NE.0 .OR. KDEC.EQ.-1)
 
135
     1    CALL TERROR(ITBL,4,'Could not find column DEC')
 
136
      CALL TBLSER (ITBL, 'EQUINOX',  KEQ,ISTAT)
 
137
      IF (ISTAT.NE.0 .OR. KEQ.EQ.-1)
 
138
     1    CALL TERROR(ITBL,5,'Could not find column EQUINOX')
 
139
C
 
140
C       Make sure table is sorted on OBJECT column:
 
141
C
 
142
      CALL TBCSRT (ITBL, 1, KOBJ, 1,  ISTAT)
 
143
      IF (ISTAT.NE.0)
 
144
     1        CALL TERROR(ITBL,8,'Could not sort OBJECT column')
 
145
C
 
146
C
 
147
C       Set up OUTPUT table file:
 
148
C
 
149
C     First, get name from tblfil local keyword:
 
150
      CALL STKRDC ('TBLFIL', 1, 1, 80,  NACTEL,TBLFIL,IUNIT,NULLS,ISTAT)
 
151
C
 
152
      CARD='Creating '//TBLFIL(:71)
 
153
      CALL TV(CARD)
 
154
      CALL TBTINI (TBLFIL, 0, 0, 1, 1,  NEWTBL, ISTAT)
 
155
C
 
156
C       create column pointers...
 
157
C
 
158
      CALL TBCINI(NEWTBL, D_C_FORMAT,32,'A32',' ','OBJECT',JOBJ,ISTAT)
 
159
      CALL TBCINI(NEWTBL, D_R4_FORMAT, 1,'R10.5',' ','RA',JRA,ISTAT)
 
160
      CALL TBCINI(NEWTBL, D_R4_FORMAT, 1,'s9.4',' ','DEC',JDEC,ISTAT)
 
161
      CALL TBCINI
 
162
     1        (NEWTBL, D_R4_FORMAT, 1,'F10.3',' ','EQUINOX',JEQ,ISTAT)
 
163
C
 
164
C     Initialize:
 
165
C
 
166
   20 ITEM=0
 
167
C
 
168
   21 NROW=NROW+1
 
169
C
 
170
C       Read row:
 
171
C
 
172
      CALL TBERDC (ITBL, NROW, KOBJ, CTEST, NULL, ISTAT)
 
173
      CALL FT_EOS (CTEST,32, OBJECT, ISTAT)
 
174
C
 
175
      IF (OBJECT.EQ.OLDOBJ .OR.OLDOBJ.EQ.' ') THEN
 
176
C         continue collecting data for this star.
 
177
          OLDOBJ=OBJECT
 
178
          ITEM=ITEM+1
 
179
          IF (ITEM.GT.MXITEM) THEN
 
180
             CALL TV('Too many data for this star.')
 
181
             CALL STETER(21,'Increase parameter MXITEM and recompile.')
 
182
          END IF
 
183
      ELSE
 
184
C         Summarize this star.
 
185
          GO TO 30
 
186
      END IF
 
187
C
 
188
      CALL TBERDR(ITBL, NROW, KRA,  RAS(ITEM),NULL, ISTAT)
 
189
      CALL TBERDR(ITBL, NROW, KDEC,  DECS(ITEM),NULL, ISTAT)
 
190
      CALL TBERDR(ITBL, NROW, KEQ,  EQUINOXS(ITEM),NULL, ISTAT)
 
191
C
 
192
      IF (NROW.EQ.NROWS) GO TO 30
 
193
      GO TO 21
 
194
C
 
195
C
 
196
C       Summarize this star:
 
197
C
 
198
   30 NROUT=NROUT+1
 
199
      ITEMS=ITEM
 
200
C
 
201
      EQUINOX=EQUINOXS(1)
 
202
      IF (ITEM.EQ.1) THEN
 
203
         RA=RAS(1)
 
204
         DEC=DECS(1)
 
205
         CARD='Only one observation of '//OLDOBJ
 
206
         CALL TVN(CARD)
 
207
         GO TO 80
 
208
      ELSE
 
209
C        extract robust estimates.
 
210
         CARD='Processing '//OLDOBJ
 
211
         CALL TVN(CARD)
 
212
      END IF
 
213
C
 
214
C       (here is where we should display scatter.)
 
215
C
 
216
      CALL SORT1(RAS,ITEMS)
 
217
      CALL SORT1(DECS,ITEMS)
 
218
C
 
219
      RA=0.5*(RAS((ITEMS+1)/2) + RAS(ITEMS/2+1))
 
220
      DEC=0.5*(DECS((ITEMS+1)/2) + DECS(ITEMS/2+1))
 
221
C
 
222
C
 
223
C        Write summary row:
 
224
C
 
225
   80 CALL TBEWRC (NEWTBL, NROUT, JOBJ, OLDOBJ,  ISTAT)
 
226
      CALL TBEWRR (NEWTBL, NROUT, JRA, RA,  ISTAT)
 
227
      CALL TBEWRR (NEWTBL, NROUT, JDEC, DEC,  ISTAT)
 
228
      CALL TBEWRR (NEWTBL, NROUT, JEQ, EQUINOX,  ISTAT)
 
229
C
 
230
      IF (NROW.EQ.NROWS) GO TO 90
 
231
C
 
232
C     Prepare for next star:
 
233
      OLDOBJ=OBJECT
 
234
      NROW=NROW-1
 
235
      GO TO 20
 
236
C
 
237
C
 
238
  90  CALL TBTCLO(ITBL,  ISTAT)
 
239
      CALL TV('   sdata.tbl closed.')
 
240
C
 
241
      CALL TBTCLO(NEWTBL,  ISTAT)
 
242
      CARD='File '//TBLFIL(1:60)//' closed'                    ! RHW 4/10/93
 
243
      CALL TV(CARD)
 
244
C
 
245
C
 
246
      CALL STSEPI
 
247
C
 
248
      END