← Back to branch summary

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

« back to all changes in this revision

Viewing changes to contrib/surfphot/libsrc/astrm.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/surfphot/libsrc/astrm.for
 
1
C @(#)astrm.for 19.1 (ES0-DMD) 02/25/03 13:30:45
 
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.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,                  
 
31
C                                         all rights reserved                   
 
32
C.IDENTIFICATION: ASTRM                                                         
 
33
C.PURPOSE:        Routine for main program star written by P.B.Stetson          
 
34
C                 performs the gaussian astrometric fitting                     
 
35
C.LANGUAGE:       F77+ESOext                                                    
 
36
C.AUTHOR:         P. Stetsan                                                    
 
37
C.ALGORITHM:      see Stetson,  P.B.,  1979  Astron. J. 84  1149                
 
38
C.INPUT/OUTPUT:   ->  KRX    input array                                        
 
39
C                 ->  JX     low  side local minimun                            
 
40
C                 ->  LX     high side local minimum                            
 
41
C                 <-> XI     input/ouput estimates of the image centre          
 
42
C                 <-  S1     error estimate of the output image centre          
 
43
C                            0. fitting procedure failed or if                  
 
44
C                            KOUNT .NE. 0  then KOUNT too small                 
 
45
C                 <-> SX     input/output estimate of the image width           
 
46
C                 <-  DO     output estimate of the image height                
 
47
C                 <-  B      output estimate of the background marginal iensity.
 
48
C                 <-> KOUNT  number of iterations required/ status              
 
49
C                            0 fitting procedure failed                         
 
50
C.COMMENTS:       If the solution fails to converge, XI will come back with     
 
51
C                 the input value.... with S1 = 0. just in case you want to try 
 
52
C                 doing the photometry with the estimated centre from SERCH.    
 
53
C.KEYWORDS :      photometry,astrometry                                         
 
54
C.VERSION:        8212?? P. Stetsan                                             
 
55
C.VERSION:        871123 R.H. Warmels  ESO Fortran Conversion                   
 
56
C.VERSION:        ?????? H. Waldhausen Vax implementation                       
 
57
C ---------------------------------------------------------------------         
 
58
      SUBROUTINE ASTRM(KRX,JX,LX,XI,S1,SX,DO,B,KOUNT)                           
 
59
C                                                                               
 
60
      IMPLICIT         NONE
 
61
      DOUBLE PRECISION KRX(1)
 
62
      INTEGER          JX
 
63
      INTEGER          LX
 
64
      REAL             XI
 
65
      REAL             S1
 
66
      REAL             SX
 
67
      REAL             DO
 
68
      REAL             B
 
69
      INTEGER          KOUNT
 
70
C
 
71
      INTEGER          I
 
72
      INTEGER          J
 
73
      INTEGER          KNT
 
74
      INTEGER          L
 
75
      REAL             D, DX, DELX
 
76
      REAL             EX
 
77
      REAL             F, FF, FOLD
 
78
      REAL             S, S2, SXSI
 
79
      REAL             XO
 
80
      REAL             Y
 
81
      DOUBLE PRECISION CTF(5)                         
 
82
      DOUBLE PRECISION CTC(5,5)                         
 
83
      DOUBLE PRECISION V(5)                         
 
84
      DOUBLE PRECISION DABS                         
 
85
      REAL             C(3)                                                     
 
86
C                                                                               
 
87
C                                                                               
 
88
C..  this part of the program uses the estimates of the image centre            
 
89
C..  and width from SERCH to derive initial estimates of the image              
 
90
C..  height and background density by linear least-squares.                     
 
91
C                                                                               
 
92
      IF ((LX-JX).GT.4) THEN                                                    
 
93
         XO     = XI                                                            
 
94
         DELX   = LX - JX + 1                                                   
 
95
         XO     = XI                                                            
 
96
         DELX   = LX - JX + 1                                                   
 
97
         S1     = 0.0                                                           
 
98
         S2     = 0.0                                                           
 
99
         DO     = 0.0                                                           
 
100
         B      = 0.0                                                           
 
101
         SXSI   = 0.5/(SX*SX)                                                   
 
102
C                                                                               
 
103
         DO 10 I = JX,LX                                                        
 
104
            DX     = I - XO                                                     
 
105
            EX     = EXP(-DX*DX*SXSI)                                           
 
106
            DO     = DO + EX*KRX(I)                                             
 
107
            B      = B + EX*EX                                                  
 
108
            S1     = S1 + KRX(I)                                                
 
109
            S2     = S2 + EX                                                    
 
110
   10    CONTINUE                                                               
 
111
C                                                                               
 
112
         F      = LX - JX + 1                                                   
 
113
         DO     = (DO-S1*S2/F)/ (B-S2*S2/F)                                     
 
114
         B      = (S1-DO*S2)/F                                                  
 
115
         IF (DO.LT.20.0) D  = 20.0                                              
 
116
         IF (SX.LT.2.0)  SX = 2.0                                               
 
117
         FOLD   = 1.0E38                                                        
 
118
C                                                                               
 
119
C..  main iteration loop                                                        
 
120
         DO 90 KNT = 1,KOUNT                           !   iteration loop       
 
121
            S1     = 0.0                                                        
 
122
            S2     = 0.0                                                        
 
123
            FF     = 0.0                                                        
 
124
            SXSI   = 1.0/ (SX*SX)                                               
 
125
            DO 30 I = 1,4                                                       
 
126
               CTF(I) = 0.0D0                                                   
 
127
               DO 20 J = I,4                                                    
 
128
                  CTC(I,J) = 0.0D0                                              
 
129
   20          CONTINUE                                                         
 
130
   30       CONTINUE                                                            
 
131
C                                                                               
 
132
C..  set normal equations                                                       
 
133
            DO 60 L = JX,LX                              !   equation loop      
 
134
               DX     = L - XO                                                  
 
135
               EX     = 0.5*DX*DX*SXSI                                          
 
136
               C(2)   = EXP(-EX)                                                
 
137
               S      = DO*C(2)                                                 
 
138
               Y      = S + B                                                   
 
139
               F      = Y - KRX(L)                                              
 
140
               C(1)   = DX*SXSI*S                                               
 
141
               EX     = C(1)*C(1)                                               
 
142
               S1     = S1 + EX*F*F                                             
 
143
               S2     = S2 + EX                                                 
 
144
               FF     = FF + F*F                                                
 
145
               C(3)   = C(1)*DX/SX                                              
 
146
               DO 50 I = 1,3                                                    
 
147
                  CTF(I) = CTF(I) + F*C(I)                                      
 
148
                  CTC(I,4) = CTC(I,4) + C(I)                                    
 
149
                  DO 40 J = I,3                                                 
 
150
                     CTC(I,J) = CTC(I,J) + C(I)*C(J)                            
 
151
   40             CONTINUE                                                      
 
152
   50          CONTINUE                                                         
 
153
               CTF(4) = CTF(4) + F  !   equation loop                           
 
154
   60       CONTINUE                                                            
 
155
C                                                                               
 
156
            CTC(4,4) = FLOAT(LX-JX+1)                                           
 
157
            DO 80 I = 2,4                                                       
 
158
               L      = I - 1                                                   
 
159
               DO 70 J = 1,L                                                    
 
160
                  CTC(I,J) = CTC(J,I)                                           
 
161
   70          CONTINUE                                                         
 
162
   80       CONTINUE                                                            
 
163
C                                                                               
 
164
C..  solve normal equations                                                     
 
165
            CALL INVRS(CTC,4)                                                   
 
166
            CALL VMUL(CTC,4,CTF,V)                                              
 
167
C                                                                               
 
168
C..  if DA is the correction to the estimate of A, then replace DA by           
 
169
C..  DA/(1+(DA/X)) = DA' prevents an ill-determined solution from               
 
170
C..  jumping around too much in subsequent iterations. DA' is less              
 
171
C..  than X, and DA' goes to DA for DA much less than X.                        
 
172
C..  If X = A/R for some R greater than one, then an intrinsically              
 
173
C..  positive parameter will never be permitted to pass through zero.           
 
174
C..  This can be helpful for parameters like image width or height.             
 
175
C                                                                               
 
176
            XO     = XO - V(1)/ (1.0+DABS(V(1)/5.0))                            
 
177
            DO     = DO - V(2)/ (1.0+2.0*DABS(V(2)/DO))                         
 
178
            SX     = SX - V(3)/ (1.0+2.0*DABS(V(3)/SX))                         
 
179
            B      = B - V(4)/ (1.0+2.0*DABS(V(4)/B))                           
 
180
C                                                                               
 
181
C..  failure                                                                    
 
182
            IF ((DO.LT.20.0) .OR. (SX.LT.1.0) .OR.                              
 
183
     +         (SX.GT.0.5*DELX)) THEN                                           
 
184
               KOUNT  = 0                                                       
 
185
               S1     = 0.0                                                     
 
186
               RETURN                                                           
 
187
            END IF                                                              
 
188
C                                                                               
 
189
C..  success                                                                    
 
190
            IF ((DABS(V(1)).LT.0.01) .AND.                                      
 
191
     +         ((FOLD-FF)/FOLD.LT.0.05)) THEN                                   
 
192
               XI     = XO                                                      
 
193
               S1     = S1*CTC(1,1)*DELX/ (S2* (DELX-4.0))                      
 
194
               RETURN                                                           
 
195
            END IF                                                              
 
196
C                                                                               
 
197
            FOLD   = FF                                                         
 
198
   90    CONTINUE                                         !   iteration loop    
 
199
C                                                                               
 
200
C..  number of iterations too small                                             
 
201
         S1     = 0.0                                                           
 
202
C                                                                               
 
203
C..  undersampled                                                               
 
204
      ELSE                                                                      
 
205
         KOUNT  = 0                                                             
 
206
         S1     = 0.0                                                           
 
207
      END IF                                                                    
 
208
C..  exit                                                                       
 
209
      RETURN                                                                    
 
210
      END