← Back to branch summary

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

« back to all changes in this revision

Viewing changes to contrib/invent/libsrc/radcor.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/radcor.for
 
1
C @(#)radcor.for        19.1 (ES0-DMD) 02/25/03 13:25:37
 
2
C===========================================================================
 
3
C Copyright (C) 1995 European Southern Observatory (ESO)
 
4
C
 
5
C This program is free software; you can redistribute it and/or 
 
6
C modify it under the terms of the GNU General Public License as 
 
7
C published by the Free Software Foundation; either version 2 of 
 
8
C the License, or (at your option) any later version.
 
9
C
 
10
C This program is distributed in the hope that it will be useful,
 
11
C but WITHOUT ANY WARRANTY; without even the implied warranty of
 
12
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
13
C GNU General Public License for more details.
 
14
C
 
15
C You should have received a copy of the GNU General Public 
 
16
C License along with this program; if not, write to the Free 
 
17
C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
 
18
C MA 02139, USA.
 
19
C
 
20
C Corresponding concerning ESO-MIDAS should be addressed as follows:
 
21
C       Internet e-mail: midas@eso.org
 
22
C       Postal address: European Southern Observatory
 
23
C                       Data Management Division 
 
24
C                       Karl-Schwarzschild-Strasse 2
 
25
C                       D 85748 Garching bei Muenchen 
 
26
C                       GERMANY
 
27
C===========================================================================
 
28
C
 
29
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
30
C.IDENTIFICATION
 
31
C  subroutine RADCOR         version 1.2        870729
 
32
C  A. Kruszewski             ESO Garching
 
33
C.PURPOSE
 
34
C  Correct different radial size parameters and isophotal magnitude
 
35
C  for influence of brightness dependent finite size of the used image.
 
36
C  A stellar object is assumed.
 
37
C  This subroutine should remove any brightness dependence of
 
38
C  corrected radial size parameters and isophotal magnitudes
 
39
C  for stellar objects.
 
40
C  It is irrelevant for galaxies.
 
41
C  It work for unsaturated objects only.
 
42
C.INPUT/OUTPUT
 
43
C  input arguments
 
44
C  PMTR        real*4 array      array holding classifiers
 
45
C  M           integer*4         number of objects
 
46
C  STPR        real*4 array      one-dimensional point spread function
 
47
C                                of stallar objects
 
48
C  TRSH        real*4            detection treshold
 
49
C  output arguments
 
50
C  PCLA        real*4 array      array holding corrected classifiers
 
51
C-----------------------------------------------------------------------
 
52
      SUBROUTINE RADCOR(PMTR,PCLA,M,STPR,TRSH)
 
53
C
 
54
      IMPLICIT NONE
 
55
      INCLUDE 'MID_REL_INCL:INVENT.INC/NOLIST'
 
56
C
 
57
      REAL     PMTR(30,MAXCNT)
 
58
      REAL     PCLA(0:13,MAXCNT)
 
59
      INTEGER  M
 
60
      REAL     STPR(MAXCNT)
 
61
      REAL     TRSH
 
62
C
 
63
      INTEGER  I, IR, IT
 
64
      INTEGER  J
 
65
      INTEGER  NPXL(0:MAXSUB)
 
66
C
 
67
      REAL     TEMP
 
68
      REAL     FI, FRCT
 
69
      REAL     SPIU(0:MAXSUB)
 
70
      REAL     DIF1, DIF2
 
71
      REAL     ALIM
 
72
      REAL     CTRS, CR, CM
 
73
      REAL     CMAG(0:100)
 
74
      REAL     CRAD(0:100),AIN0(0:MAXSUB)
 
75
      REAL     AIN1(0:MAXSUB),AIN2(0:MAXSUB),AINK(0:MAXSUB)
 
76
C
 
77
      DATA NPXL/1,8,12,16,32,28,40,40,48,68,56,72,68,88,88,84,112,112,
 
78
     +     112,116,112,144,140,144,144,168,164,160,184,172,200,192,188,
 
79
     +     208,224,224,228,224,248,236,264,248,264,276,264,288,276,304,
 
80
     +     304,312,316/
 
81
C
 
82
C     Transforms the differential logarithmic
 
83
C *** point spread function STPR into profile
 
84
C     in intensity units SPIU.
 
85
C
 
86
      SPIU(0) = 1.0
 
87
      DO 10 I = 1,50
 
88
          SPIU(I) = SPIU(I-1)* (10.0** (-STPR(I)))
 
89
   10 CONTINUE
 
90
C
 
91
C *** Calculates integrated quantities.
 
92
C
 
93
      AIN0(0) = 1.0
 
94
      AIN1(0) = 0.333
 
95
      AIN2(0) = 0.1
 
96
      AINK(0) = 10.0
 
97
      DO 20 I = 1,50
 
98
          TEMP   = SPIU(I)*NPXL(I)
 
99
          AIN0(I) = AIN0(I-1) + TEMP
 
100
          AIN1(I) = AIN1(I-1) + TEMP*I
 
101
          AIN2(I) = AIN2(I-1) + TEMP*I*I
 
102
          AINK(I) = AINK(I-1) + TEMP/ (I*I)
 
103
   20 CONTINUE
 
104
      DO 30 I = 0,50
 
105
          AIN1(I) = AIN1(I)/AIN0(I)
 
106
          AIN2(I) = AIN2(I)/AIN0(I)
 
107
          AINK(I) = AINK(I)/AIN0(I)
 
108
   30 CONTINUE
 
109
      DO 40 I = 0,50
 
110
          AIN1(I) = AIN1(25)/AIN1(I)
 
111
          AIN2(I) = AIN2(25)/AIN2(I)
 
112
          AIN2(I) = AIN2(I)**0.5
 
113
          AINK(I) = AINK(25)/AINK(I)
 
114
          AINK(I) = AINK(I)**(-0.5)
 
115
   40 CONTINUE
 
116
C
 
117
C *** Correct following parameters: "RADIUS1", "RADIUS2","KRON RADIUS".
 
118
C
 
119
      DO 50 I = 1,M
 
120
C
 
121
C *** Store value of the limiting radius in variable TEMP.
 
122
          TEMP   = PMTR(15,I)
 
123
          IR     = INT(TEMP)
 
124
          FRCT   = TEMP - IR
 
125
C
 
126
C *** Apply the corrections.
 
127
C
 
128
          PCLA(4,I)  = PMTR(6,I)*  
 
129
     2                 (AIN1(IR)* (1.0-FRCT)+AIN1(IR+1)*FRCT)
 
130
          PCLA(5,I)  = PMTR(7,I)*  
 
131
     2                 (AIN2(IR)* (1.0-FRCT)+AIN2(IR+1)*FRCT)
 
132
          PCLA(12,I) = PMTR(18,I)* 
 
133
     2                 (AINK(IR)* (1.0-FRCT)+AINK(IR+1)*FRCT)
 
134
   50 CONTINUE
 
135
C
 
136
C     Corrects parameters: "ISOPHOTAL MAGNITUDE", "LIMITING RADIUS".
 
137
C     it takes into account that the stellar profile is cut by
 
138
C     limiting treshold in a different way depending on stellar
 
139
C *** brightness. Average of 9 central pixels is taken as a 
 
140
C     measure of stellar brightness. Stellar profile is cut 
 
141
C     horizontally by 100 eqidistant levels. Corrections are
 
142
C     calculated for each level and stored in arrays CMAG and 
 
143
C     CRAD. Linear interpolation is used.
 
144
C
 
145
      DO 80 I = 1,100
 
146
          J      = 0
 
147
          CMAG(I) = 0.0
 
148
          CRAD(I) = 0.0
 
149
          FI     = FLOAT(I)
 
150
          ALIM   = (100.0-FI)*0.01
 
151
   60     CONTINUE
 
152
              IF ( .NOT. (SPIU(J).GT.ALIM) ) GO TO 70
 
153
              CMAG(I) = CMAG(I) + SPIU(J)*NPXL(J)
 
154
              J = J + 1
 
155
              IF (J.GT.50) THEN
 
156
                 J = 50
 
157
              ELSE
 
158
                GO TO 60
 
159
              ENDIF
 
160
   70     CONTINUE
 
161
          DIF1   = SPIU(J-1) - ALIM
 
162
          DIF2   = SPIU(J-1) - SPIU(J)
 
163
          IF (DIF1.GT.0.0001 .AND. DIF2.GT.0.0001) THEN
 
164
              FRCT   = DIF1/DIF2
 
165
          ELSE
 
166
              FRCT   = 0.0
 
167
          END IF
 
168
 
 
169
          IF (FRCT.GT.1.0) THEN
 
170
              FRCT   = 1.0
 
171
          END IF
 
172
 
 
173
          CMAG(I) = CMAG(I) + FRCT*SPIU(J)*NPXL(J)
 
174
          CRAD(I) = J - 1 + FRCT
 
175
   80 CONTINUE
 
176
C
 
177
      DO 90 I = 1,100
 
178
          CMAG(I) = 2.5*LOG10(CMAG(I)/CMAG(100))
 
179
          CRAD(I) = CRAD(99) - CRAD(I)
 
180
   90 CONTINUE
 
181
C
 
182
      CMAG(0) = 3.0*CMAG(1) - 2.0*CMAG(2)
 
183
      CRAD(0) = 3.0*CRAD(1) - 2.0*CRAD(2)
 
184
      CTRS   = (SPIU(0)+8.0*SPIU(1))* (TRSH/9)
 
185
      DO 100 I = 1,M
 
186
          IF (PMTR(2,I).GT.TRSH) THEN
 
187
              TEMP   = (1.0-CTRS/PMTR(2,I))*100.0
 
188
              IT     = INT(TEMP)
 
189
              FRCT   = TEMP - IT
 
190
              CM     = CMAG(IT)* (1.0-FRCT) + CMAG(IT+1)*FRCT
 
191
              IF (IT.EQ.99) THEN
 
192
                  CR     = 0.0
 
193
 
 
194
              ELSE
 
195
                  CR     = CRAD(IT)* (1.0-FRCT) + CRAD(IT+1)*FRCT
 
196
              END IF
 
197
 
 
198
          ELSE
 
199
              CM     = 0.0
 
200
              CR     = 0.0
 
201
          END IF
 
202
 
 
203
          PCLA(3,I) = PMTR(5,I) + CM
 
204
          PCLA(9,I) = PMTR(15,I) + CR
 
205
  100 CONTINUE
 
206
C
 
207
      RETURN
 
208
      END