← Back to branch summary

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

« back to all changes in this revision

Viewing changes to contrib/invent/libsrc/cmpsub.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/invent/libsrc/cmpsub.for
 
1
C===========================================================================
 
2
C Copyright (C) 1995 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 Massachusetss Ave, Cambridge, 
 
17
C MA 02139, USA.
 
18
C
 
19
C Corresponding 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
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
29
C.IDENT:   subroutine CMPSUB          version 1.1          870315
 
30
C          A. Kruszewski              Obs. de Geneve
 
31
C          undersampled case          version 2.0          880826
 
32
C          A. Kruszewski              Warsaw U. Obs.
 
33
C.PURPOSE: Subtracts stellar components in a subframe copied from image frame.
 
34
C.INPUT/OUTPUT
 
35
C  Input parameters
 
36
C
 
37
C  A         real*4 array        input image array
 
38
C  JAPY      integer*4 array     pointers to imege lines
 
39
C  IBUF      integer*4           limits of image buffer
 
40
C  I         integer*4           x-coordinate of object
 
41
C  J         integer*4           y-coordinate of object
 
42
C  L0        integer*4           offset to catalog objects
 
43
C  L1        integer*4           last object in catalog
 
44
C  LW        integer*4           last object written to file
 
45
C  L         integer*4           identification number of object
 
46
C  M         integer*4           actuall number of objects
 
47
C  LSTP      integer*4 array     manage linked lists of pointers
 
48
C  NREG      integer*4           number of regional linked lists
 
49
C  LHED      integer*4           half edge size of subarray
 
50
C  TRSH      real*4              limiting threshold
 
51
C  APSF      real*4 array        one-dimensional point spread function
 
52
C  FPSF      real*4 array        two-dimensional point spread function
 
53
C  IPSF      integer*4 array     pointers to two-dimensional p.s.f.
 
54
C  LPXL      integer*4           number of 2-d pixels
 
55
C  LSBP      integer*4           number of subpixels
 
56
C  NCAT      integer*4 array     integer parameters array
 
57
C  PMTR      real*4 array        real parameters array
 
58
C  PRCT      real*4 array        catalog of one-dimensional profiles
 
59
C  NP        integer*4           dimension of catalog buffer
 
60
C  NPAS      integer*4           actuall iteration number
 
61
C  HCUT      real*4              saturation level
 
62
C  Output parameters
 
63
C  AS        real*4 array        data subarray after subtraction
 
64
C  JAPYS     integer*4 array     pointers to subarray lines
 
65
C  IBUFS     integer*4 array     limits of data in subarray
 
66
C----------------------------------------------------------------------
 
67
      SUBROUTINE CMPSUB (A, JAPY, IBUF, I, J, L0, L1, LW,
 
68
     &            L, M, NREG, LSTP, LHED, TRSH, APSF, FPSF,
 
69
     &            IPSF, LPXL, LSBP, NCAT, PMTR, PRCT, NPAS,
 
70
     &            AS, JAPYS, IBUFS, HCUT, MAS )
 
71
C
 
72
      IMPLICIT  NONE
 
73
      INCLUDE  'MID_REL_INCL:INVENT.INC/NOLIST'
 
74
C
 
75
      INTEGER   NP , NREG
 
76
      INTEGER   JAPY(1) , IBUF(4), NCAT(NIPAR,MAXCNT)
 
77
      INTEGER   JAPYS(2*MAXSUB+1), IBUFS(4)
 
78
      INTEGER   LSTP(0:4, 0:NREG)
 
79
      INTEGER   NCT(NIPAR)
 
80
      INTEGER   MAS((-MAXSUB):MAXSUB,(-MAXSUB):MAXSUB)
 
81
      INTEGER   I, IPSF(1), J, L, L0, L1, LHED
 
82
      INTEGER   LPXL, LSBP, LW
 
83
      INTEGER   M
 
84
      INTEGER   NPAS
 
85
      INTEGER   MDIS
 
86
      INTEGER   MPXL
 
87
      INTEGER   MSBP
 
88
      INTEGER   IL1, IL2
 
89
      INTEGER   JL1, JL2
 
90
      INTEGER   MXS
 
91
      INTEGER   ITMP
 
92
      INTEGER   IS, IS1, IS2
 
93
      INTEGER   JOFF, JOFFS
 
94
      INTEGER   JOF, JOFS
 
95
      INTEGER   IL, IN
 
96
      INTEGER   K, KK, KK2
 
97
      INTEGER   LK, LN, LZ, LD
 
98
      INTEGER   JN, JL
 
99
      INTEGER   IDIF
 
100
      INTEGER   JDIF
 
101
      INTEGER   IOPC
 
102
      INTEGER   IINT
 
103
C
 
104
      REAL      AS((2*MAXSUB+1)**2)
 
105
      REAL      A(1), APSF(0:MAXSUB)
 
106
      REAL      CINT
 
107
      REAL      CINT1
 
108
      REAL      FPSF(1)
 
109
      REAL      HCUT
 
110
      REAL      PRCT( 0:MAXSUB, MAXCNT), PMTR(NRPAR,MAXCNT)
 
111
      REAL      PMT(NRPAR), PCT(0:MAXSUB)
 
112
      REAL      SCALE
 
113
      REAL      TRSH
 
114
C
 
115
C     LOGICAL   CENTER
 
116
      LOGICAL   UNDER
 
117
      LOGICAL   DONE
 
118
C
 
119
      IINT = 0
 
120
      NP = MAXCNT
 
121
      MDIS = MAXSUB + LHED
 
122
      IF ( LPXL .GT. 0 .OR. LSBP .GT. 0 ) THEN
 
123
            UNDER = .TRUE.
 
124
            MPXL = 2 * LPXL + 1
 
125
            MSBP = 2 * LSBP + 1
 
126
      ELSE
 
127
            UNDER = .FALSE.
 
128
      ENDIF
 
129
      LHED = MIN( MAXSUB , LHED )
 
130
      SCALE = 9.0 / ( 1.0 + 8.0*APSF(1) )
 
131
      CINT = 0.0
 
132
C
 
133
C ******      Prepare subarrays.
 
134
C
 
135
      IL1 = MAX( I-LHED , IBUF(1) )
 
136
      IL2 = MIN( I+LHED , IBUF(3) )
 
137
      JL1 = MAX( J-LHED , IBUF(2) )
 
138
      JL2 = MIN( J+LHED , IBUF(4) )
 
139
      IBUFS(1) = IL1 - I
 
140
      IBUFS(2) = JL1 - J
 
141
      IBUFS(3) = IL2 - I
 
142
      IBUFS(4) = JL2 - J
 
143
      MXS = 2*MAXSUB + 1
 
144
      IS1 = IBUFS(2)
 
145
      IS2 = IBUFS(4)
 
146
      DO 5 IS = IS1 , IS2
 
147
          ITMP = IS - IS1 + 1
 
148
          JAPYS(ITMP) = (ITMP-1) * MXS + MAXSUB + 1
 
149
    5 CONTINUE
 
150
      JOFF = IBUF(2) - 1
 
151
      JOFFS = IBUFS(2) - 1
 
152
      DO 10 JL = JL1 , JL2
 
153
          JOF = JAPY(JL-JOFF)
 
154
          JOFS = JAPYS(JL-J-JOFFS)
 
155
          DO 20 IL = IL1 , IL2
 
156
              AS(JOFS+IL-I) = A(JOF+IL)
 
157
   20     CONTINUE
 
158
   10 CONTINUE
 
159
C
 
160
C *** Start subtracting contributions from close components.
 
161
C
 
162
      K = 0
 
163
   31 CONTINUE
 
164
          CALL GETLST( L , L0 , L1 , MDIS , NREG ,
 
165
     &                 LSTP , NCAT , PMTR , PRCT , K ,
 
166
     &                 NCT , PMT , PCT , DONE )
 
167
          IF ( .NOT. DONE ) GOTO 30
 
168
C
 
169
C ***     Ignore the L-th object.
 
170
C
 
171
          IF ( K .EQ. L ) GOTO 32
 
172
C
 
173
C ***     Look for close neighbour.
 
174
C ***     Check first if K-th object is in buffer.
 
175
C
 
176
c            IF ( K .LE. L0 .OR. K .GT. L1 ) THEN
 
177
c                  READ( ISF , REC=K )  NCT , PMT , PCT
 
178
c            ELSE
 
179
          LK = K - L0
 
180
          DO 60 KK = 1 , NIPAR
 
181
              NCT(KK) = NCAT(KK,LK)
 
182
   60     CONTINUE
 
183
          DO 70 KK = 1 , NRPAR
 
184
              PMT(KK) = PMTR(KK,LK)
 
185
   70     CONTINUE
 
186
          KK2 = MIN( NCT(5) , MAXSUB )
 
187
          DO 80 KK = 0 , KK2
 
188
                PCT(KK) = PRCT(KK,LK)
 
189
   80     CONTINUE
 
190
c         ENDIF
 
191
          IF ( PMT(2) .LE. TRSH ) THEN
 
192
C
 
193
C ***         Too faint objects are not considered.
 
194
C
 
195
              GOTO 32
 
196
          ENDIF
 
197
          IN = NCT(1)
 
198
          JN = NCT(2)
 
199
          LN = NCT(5)
 
200
          IDIF = I-IN
 
201
          JDIF = J-JN
 
202
          LD = LN + LHED
 
203
          IF ( ABS(IDIF) .LE. LD .AND. ABS(JDIF) .LE. LD ) THEN
 
204
C
 
205
C ***         Found one!
 
206
C
 
207
              IOPC = -1
 
208
              CALL STARSA( IOPC , AS , JAPYS , IBUFS , I ,
 
209
     &                     J , LPXL , LSBP , NCT , PMT ,
 
210
     &                     PCT , APSF , FPSF , NPAS , SCALE ,
 
211
     &                     CINT1 )
 
212
              IF ( CINT1 .GT. CINT ) THEN
 
213
                  CINT = CINT1
 
214
                  IINT = K
 
215
              ENDIF
 
216
          ENDIF
 
217
   32     CONTINUE
 
218
          IF ( DONE ) GOTO 31
 
219
   30 CONTINUE
 
220
C
 
221
C *** Write down the number of the brightest stellar component
 
222
C *** if it contribute more than 10 % to the central pixel.
 
223
C
 
224
      LZ = L - L0
 
225
      IF ( CINT .GT. 0.1*PMTR(2,LZ) ) THEN
 
226
          NCAT(10,LZ) = IINT
 
227
      ELSE
 
228
          NCAT(10,LZ) = 0
 
229
      ENDIF
 
230
C
 
231
C *** Restore values of saturated points.
 
232
C
 
233
      DO 40 JL = JL1 , JL2
 
234
          JOFS = JAPYS(JL-J-JOFFS)
 
235
          DO 50 IL = IL1 ,IL2
 
236
              IF ( MAS(IL-I,JL-J) .EQ. -1 ) THEN
 
237
                  AS(JOFS+IL-I) = HCUT
 
238
              ENDIF
 
239
   50     CONTINUE
 
240
   40 CONTINUE
 
241
C
 
242
      RETURN
 
243
C
 
244
      END