← Back to branch summary

~ubuntu-branches/ubuntu/wily/eso-midas/wily-proposed

« back to all changes in this revision

Viewing changes to prim/plot/libsrc/futils.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:
prim/plot/libsrc/futils.for
 
1
C===========================================================================
 
2
C Copyright (C) 1995-2008 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
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
29
C.IDENTIFICATION: FUTILS
 
30
C.PURPOSE:        subroutines used by centerrow.for
 
31
C.LANGUAGE:       F77+ESOext
 
32
C.AUTHOR:         J.D.Ponz
 
33
C.KEYWORDS:       Line, centre
 
34
C.ALGORITHM:      Line center is found by:
 
35
C                 maximum/minimum - position of the maximum/minimum value
 
36
C                 gaussian        - center of the fitted gaussian
 
37
C                 gravity         - gravity center of the 2 highest pixels
 
38
C                                   with respect to the third one
 
39
C 
 
40
C.VERSION: 840328  JDP  Creation   
 
41
C.VERSION: 071016  KB   copied from ../src/centerrow.for
 
42
C 
 
43
C 080724        last modif
 
44
C-----------------------------------------------------------------
 
45
C 
 
46
      SUBROUTINE PLFIND(X,XSTR,XSTP,IPL1,IPL2,IMODE,IMETH,
 
47
     +                XCENTR,PEAK,IFAIL,W1,W2,ACOE,Y1,Y2)
 
48
C++++
 
49
C. PURPOSE: Find center
 
50
C---
 
51
      IMPLICIT NONE
 
52
      INTEGER          IPL1,IPL2,IMODE,IMETH
 
53
      INTEGER          IFAIL,NPX
 
54
 
 
55
      REAL             X(*)
 
56
      DOUBLE PRECISION XSTR,XSTP
 
57
      DOUBLE PRECISION XCENTR,PEAK,CHICHK,XSP,XGO
 
58
      DOUBLE PRECISION W1(*),W2(*),ACOE(*),Y1,Y2
 
59
C
 
60
      CHICHK = 0.005
 
61
      XSP    = XSTP
 
62
      XGO    = XSTR + (IPL1-1)*XSP
 
63
      NPX    = IPL2 - IPL1 + 1
 
64
      Y1     = DBLE(X(IPL1))
 
65
      Y2     = DBLE(X(IPL2))
 
66
      IF (IMETH.LT.0) THEN                                    ! gaussian center
 
67
         CALL SGAUS(X(IPL1),W1,W2,IMODE,NPX,IFAIL,XGO,XSP,
 
68
     +              XCENTR,CHICHK,0,PEAK,ACOE)
 
69
      ELSE IF (IMETH.EQ.0) THEN
 
70
         CALL GRAVT(X(IPL1),NPX,IMODE,IFAIL,XGO,XSP,XCENTR,PEAK)
 
71
      ELSE
 
72
         CALL CNTRH(X(IPL1),NPX,IMODE,IFAIL,XGO,XSP,XCENTR,PEAK)
 
73
      ENDIF
 
74
C 
 
75
      RETURN
 
76
      END
 
77
 
 
78
C 
 
79
      SUBROUTINE GRAVT(WINDOW,NPIX,IWID,IWERR,XGO,DXSTP,XCNTR,AM)       
 
80
C                                                                              
 
81
C line centering by finding the center of gravity of 2 highest                 
 
82
C points of 3 points                                                           
 
83
C the 2 highest points have values relatives to the third                      
 
84
C only applicable in emission mode                                             
 
85
C PARAMETERS :                                                                 
 
86
C WINDOW buffer with expected line center in middle                            
 
87
C NPIX number of samples in WINDOW                                             
 
88
C IWID code 0 absorption, 1 emission                                           
 
89
C IWERR  error return, 0 ok, 1 error                                           
 
90
C XGO xcoordinate first element in WINDOW                                      
 
91
C DXSTP step in x                                                              
 
92
C XCNTR return x center calculated                                             
 
93
C AM intensity center pixel                                                    
 
94
C                                                                              
 
95
      IMPLICIT NONE
 
96
C
 
97
      INTEGER  NPIX, IWID, IWERR, I, K
 
98
 
 
99
      REAL             WINDOW(*)
 
100
      DOUBLE PRECISION AM, AL, AH, DR
 
101
      DOUBLE PRECISION XGO, DXSTP, XCNTR, A, B, XSHIFT
 
102
C
 
103
      IF (IWID.NE.1) GO TO 20                                                  
 
104
C                                                                              
 
105
C find maximum                                                                 
 
106
C                                                                              
 
107
      AM     = WINDOW(1)                                                       
 
108
      I      = 1                                                               
 
109
      DO 10 K = 2,NPIX                                                         
 
110
          IF (WINDOW(K).GT.AM) THEN                                            
 
111
              AM     = WINDOW(K)                                               
 
112
              I      = K                                                       
 
113
          END IF                                                               
 
114
   10 CONTINUE                                                                 
 
115
C                                                                              
 
116
C check not boundary                                                           
 
117
C                                                                              
 
118
      IF (I.EQ.1 .OR. I.EQ.NPIX) GO TO 20                                      
 
119
      XCNTR  = XGO + (I-1)*DXSTP                                               
 
120
C                                                                              
 
121
C find lowest of the flanking pixels                                           
 
122
C                                                                              
 
123
      AL     = WINDOW(I-1)                                                     
 
124
      AH     = WINDOW(I+1)                                                     
 
125
      DR     = 1.                                                              
 
126
      IF (AL.GE.AH) THEN                                                       
 
127
          AL     = WINDOW(I+1)                                                 
 
128
          AH     = WINDOW(I-1)                                                 
 
129
          DR     = -1.                                                         
 
130
      END IF                                                                   
 
131
      AM     = WINDOW(I)                                                       
 
132
      A      = AM - AL                                                         
 
133
      B      = AH - AL                                                         
 
134
      XSHIFT = (B/ (A+B))*DXSTP                                                
 
135
      XCNTR  = XCNTR + XSHIFT*DR                                               
 
136
      IWERR  = 0                                                               
 
137
      RETURN                                                                   
 
138
C                                                                              
 
139
C error return                                                                 
 
140
C                                                                              
 
141
   20 IWERR  = 1                                                               
 
142
      RETURN                                                                   
 
143
                                                                               
 
144
      END
 
145
 
 
146
 
 
147
      SUBROUTINE CNTRH(WINDOW,NPIX,IWID,IWERR,XGO,DXSTP,XCNTR,AM)              
 
148
C                                                                              
 
149
C center of line found by simple maximum within a window                       
 
150
C                                                                              
 
151
C PARAMETERS :                                                                 
 
152
C WINDOW buffer with expected line center in middle                            
 
153
C NPIX number of samples in WINDOW                                             
 
154
C IWID code 0 absorption, 1 emission                                           
 
155
C IWERR  error return, 0 ok, 1 error                                           
 
156
C XGO xcoordinate first element in WINDOW                                      
 
157
C DXSTP step in x                                                              
 
158
C XCNTR return x center calculated                                             
 
159
C AM intensity center pixel                                                    
 
160
C                                                                              
 
161
      IMPLICIT NONE
 
162
C
 
163
      INTEGER           NPIX, IWID, IWERR
 
164
      REAL              WINDOW(*)
 
165
      DOUBLE PRECISION  AM, XGO, DXSTP, XCNTR
 
166
C
 
167
      INTEGER IPN, I
 
168
C
 
169
      IWERR  = 0                                                               
 
170
      AM     = WINDOW(1)                                                       
 
171
      IPN    = 1                                                               
 
172
      IF (IWID.NE.1) THEN                                                      
 
173
C                                                                              
 
174
C find minimum                                                                 
 
175
C                                                                              
 
176
          DO 10 I = 2,NPIX                                                     
 
177
              IF (WINDOW(I).LT.AM) THEN                                        
 
178
                  AM     = WINDOW(I)                                           
 
179
                  IPN    = I                                                   
 
180
              END IF                                                           
 
181
   10     CONTINUE                                                             
 
182
      ELSE                                                                     
 
183
C                                                                              
 
184
C find maximum                                                                 
 
185
C                                                                              
 
186
          DO 20 I = 2,NPIX                                                     
 
187
              IF (WINDOW(I).GT.AM) THEN                                        
 
188
                  AM     = WINDOW(I)                                           
 
189
                  IPN    = I                                                   
 
190
              END IF                                                           
 
191
   20     CONTINUE                                                             
 
192
      END IF                                                                   
 
193
C                                                                              
 
194
C check result                                                                 
 
195
C                                                                              
 
196
      IF (IPN.EQ.1 .OR. IPN.EQ.NPIX) GO TO 30                                  
 
197
      XCNTR  = XGO + (IPN-1)*DXSTP                                             
 
198
      RETURN                                                                   
 
199
C                                                                              
 
200
C error return                                                                 
 
201
C                                                                              
 
202
   30 IWERR  = 1                                                               
 
203
      RETURN                                                                   
 
204
      END                                                                      
 
205
 
 
206
      SUBROUTINE SGAUS(WINDOW,XBUF,YFIT,IWID,NPIX,IWERR,XGO,DXSTP,
 
207
     2                 XCNTR,CHICHK,IOUTPT,AM,A)                              
 
208
C                                                                              
 
209
C line center finding with gaussian fit                                        
 
210
C                                                                              
 
211
C PARAMETERS :                                                                 
 
212
C WINDOW buffer with expected line center in middle                            
 
213
C XBUF work buffer                                                             
 
214
C YFIT work buffer                                                             
 
215
C IWID code 0 absorption, 1 emission                                           
 
216
C NPIX number of samples in WINDOW                                             
 
217
C IWERR  error return, 0 ok, 1 error                                           
 
218
C XGO xcoordinate first element in WINDOW                                      
 
219
C DXSTP step in x                                                              
 
220
C XCNTR return x center calculated                                             
 
221
C CHICHK check chisquare                                                       
 
222
C IOUTPT display intermediate results 0 no, 1 yes                              
 
223
C AM intensity center pixel                                                    
 
224
C A array(4) coefficients                                                      
 
225
C                                                                              
 
226
      INTEGER          I,K1,K2,NPIX,IWID,IWERR,IOUTPT
 
227
      INTEGER          ITCNT,ITCHK,JER,IBRK
 
228
 
 
229
      REAL             WINDOW(*)
 
230
      DOUBLE PRECISION XBUF(*),YFIT(*),A(*)          
 
231
      DOUBLE PRECISION DELTAA(4),SIGMAA(4),XCNTR,XGO,CHICHK
 
232
      DOUBLE PRECISION DXSTP,CHISQR,FLAMDA,OLDCH,DIF,HALF,WIDTH
 
233
      DOUBLE PRECISION AM,X,STOP,TOPM
 
234
C                                                                              
 
235
C buffer with x values                                                         
 
236
C                                                                              
 
237
      X      = XGO                                                             
 
238
      DO 10 I = 1,NPIX                                                         
 
239
          XBUF(I) = XGO + (I-1)*DXSTP                                          
 
240
   10 CONTINUE                                                                 
 
241
C                                                                              
 
242
C find initial guess                                                           
 
243
C                                                                              
 
244
      CALL CNTRH(WINDOW,NPIX,IWID,IWERR,XGO,DXSTP,A(2),AM)                     
 
245
      IF (IWERR.NE.0) GO TO 120                                                
 
246
      A(4)   = (WINDOW(1)+WINDOW(NPIX))/2.                                     
 
247
      A(1)   = AM - A(4)                                                       
 
248
C                                                                              
 
249
C calculate half width                                                         
 
250
C                                                                              
 
251
      HALF   = A(1)/2. + A(4)                                                  
 
252
C                                                                              
 
253
C calculation of FWHM depends on sign of the line                              
 
254
C                                                                              
 
255
      IF (IWID.NE.1) THEN                                                      
 
256
          DO 20 I = 1,NPIX                                                     
 
257
              IF (WINDOW(I).LT.HALF) GO TO 30                                  
 
258
   20     CONTINUE                                                             
 
259
   30     K1     = I                                                           
 
260
          DO 40 I = K1,NPIX                                                    
 
261
              IF (WINDOW(I).GT.HALF) GO TO 50                                  
 
262
   40     CONTINUE                                                             
 
263
   50     K2     = I                                                           
 
264
                                                                               
 
265
      ELSE                                                                     
 
266
          DO 60 I = 1,NPIX                                                     
 
267
              IF (WINDOW(I).GT.HALF) GO TO 70                                  
 
268
   60     CONTINUE                                                             
 
269
   70     K1     = I                                                           
 
270
          DO 80 I = K1,NPIX                                                    
 
271
              IF (WINDOW(I).LT.HALF) GO TO 90                                  
 
272
   80     CONTINUE                                                             
 
273
   90     K2     = I                                                           
 
274
      END IF                                                                   
 
275
                                                                               
 
276
      WIDTH  = ABS(FLOAT(K2-K1)*DXSTP)                                         
 
277
      A(3)   = WIDTH/2.35482                !sigma = fwhm/2.35482
 
278
      OLDCH  = 9.E16                                                           
 
279
C                                                                              
 
280
C call fitting routine                                                         
 
281
C iterate 50 times                                                             
 
282
C                                                                              
 
283
      ITCNT  = 0                                                               
 
284
      ITCHK  = 50                                                              
 
285
  100 FLAMDA = 0.001                                                           
 
286
      IBRK   = 0                                                               
 
287
C                                                                              
 
288
      CALL CURFI(XBUF,WINDOW,SIGMAA,NPIX,4,0,A,DELTAA,FLAMDA,YFIT,             
 
289
     +           CHISQR,JER)                                                   
 
290
      IF (JER.NE.0) GO TO 120                                                  
 
291
C                                                                              
 
292
C calc and display intermediate results                                        
 
293
C                                                                              
 
294
      IF (IBRK.EQ.-1) GO TO 110                                                
 
295
C                                                                              
 
296
C see if CHISQR has changed significantly                                      
 
297
C                                                                              
 
298
      DIF    = OLDCH - CHISQR                                                  
 
299
      STOP   = DIF/CHISQR                                                      
 
300
      OLDCH  = CHISQR                                                          
 
301
      ITCNT  = ITCNT + 1                                                       
 
302
      IF (ITCNT.GT.ITCHK) GO TO 120                                            
 
303
      IF (STOP.GT.CHICHK) GO TO 100                                            
 
304
  110 XCNTR  = A(2)                                                            
 
305
C                                                                              
 
306
C last check that we are within the window                                     
 
307
C                                                                              
 
308
      TOPM   = XBUF(NPIX)                                                      
 
309
      IF (DXSTP.LT.0.) THEN                                                    
 
310
          IF ((XCNTR.GT.XGO) .OR. (XCNTR.LT.TOPM)) GO TO 120                   
 
311
                                                                               
 
312
      ELSE                                                                     
 
313
          IF ((XCNTR.LT.XGO) .OR. (XCNTR.GT.TOPM)) GO TO 120                   
 
314
      END IF                                                                   
 
315
                                                                               
 
316
      IWERR  = 0                                                               
 
317
      A(3)   = A(3)*2.35482                  !fwhm = sigma * 2.35482
 
318
      RETURN                                                                   
 
319
C                                                                              
 
320
C error return                                                                 
 
321
C                                                                              
 
322
  120 IWERR  = 1                                                               
 
323
      RETURN                                                                   
 
324
                                                                               
 
325
      END                                                                      
 
326
 
 
327
      SUBROUTINE CURFI(X,Y,SIGMAY,NPTS,NTERMS,MODE,A,DELTAA,FLAMDA,YFIT,       
 
328
     +                 CHISQR,IERR)                                            
 
329
C                                                                              
 
330
C least squares fit to a non-linear function                                   
 
331
C                                                                              
 
332
C PARAMETERS:                                                                  
 
333
C X array of data points ind. var.                                             
 
334
C Y array of data points dep. var.                                             
 
335
C SIGMAY array of std dev for Y data points                                    
 
336
C NPTS number of pairs of data points                                          
 
337
C MODE determines method of weighting least-squares fit                        
 
338
C  +1 (instrumental) weight = 1./sigmay(i)**2                                  
 
339
C   0 (no weighting) weight = 1.                                               
 
340
C  -1 (statistical ) weight = 1./y(i)                                          
 
341
C A array of parameters                                                        
 
342
C DELTAA array of increments for parameters A                                  
 
343
C SIGMAA array of standard deviations for parameters A                         
 
344
C FLAMDA proportion of gradient search included                                
 
345
C YFIT array of calculated values of Y                                         
 
346
C CHISQR reduced chi square for fit                                            
 
347
C IERR error return 0 ok, 1 not converging                                     
 
348
C                                                                              
 
349
C dimension statement valid for NTERMS up to 10                                
 
350
C set FLAMDA to 0.001 at the beginning of the search                           
 
351
C     
 
352
      INTEGER I,J,K,NPTS,NFREE ,NTERMS,IERR  
 
353
      INTEGER MODE ,ICNT, ISTAT                               
 
354
 
 
355
      REAL             Y(*)
 
356
      DOUBLE PRECISION X(*), SIGMAY(*), DELTAA(*), YFIT(*) 
 
357
      DOUBLE PRECISION A(*), DET, FLAMDA, CHISF, FUNCT 
 
358
      DOUBLE PRECISION WEIGHT(2048), ALPHA(10,10), BETA(10) 
 
359
      DOUBLE PRECISION DERIV(10), B(10), CHISQR, CHISQ1 
 
360
      DOUBLE PRECISION ARRAY(10,10)
 
361
      EXTERNAL FUNCT
 
362
C                                                                              
 
363
      ICNT   = 0                                                               
 
364
      IERR   = 1                                                               
 
365
      NFREE  = NPTS - NTERMS                                                   
 
366
      IF (NFREE) 20,20,30                                                      
 
367
   20 CHISQR = 0.                                                              
 
368
      RETURN                                                                   
 
369
C                                                                              
 
370
C evaluate weights                                                             
 
371
C                                                                              
 
372
   30 CONTINUE                                                                 
 
373
      IERR   = 0                                                               
 
374
      DO 100 I = 1,NPTS                                                        
 
375
          IF (MODE) 40,70,80                                                   
 
376
   40     IF (Y(I)) 60,70,50                                                   
 
377
   50     WEIGHT(I) = 1./Y(I)                                                  
 
378
          GO TO 90                                                             
 
379
                                                                               
 
380
   60     WEIGHT(I) = 1./ (-Y(I))                                              
 
381
          GO TO 90                                                             
 
382
                                                                               
 
383
   70     WEIGHT(I) = 1.                                                       
 
384
          GO TO 90                                                             
 
385
                                                                               
 
386
   80     WEIGHT(I) = 1./SIGMAY(I)**2                                          
 
387
   90     CONTINUE                                                             
 
388
  100 CONTINUE                                                                 
 
389
C                                                                              
 
390
C evaluate alpha and beta matrices                                             
 
391
C                                                                              
 
392
      DO 120 J = 1,NTERMS                                                      
 
393
          BETA(J) = 0.                                                         
 
394
          DO 110 K = 1,J
 
395
              ALPHA(J,K) = 0.                                                  
 
396
  110     CONTINUE                                                             
 
397
  120 CONTINUE                                                                 
 
398
      DO 150 I = 1,NPTS                                                        
 
399
          CALL FDERI(X,I,A,DELTAA,NTERMS,DERIV)                                
 
400
          DO 140 J = 1,NTERMS                                                  
 
401
              BETA(J) = BETA(J) + WEIGHT(I)* (Y(I)-FUNCT(X,I,A))*              
 
402
     +                  DERIV(J)                                               
 
403
              DO 130 K = 1,J                                                   
 
404
                  ALPHA(J,K) = ALPHA(J,K) + WEIGHT(I)*DERIV(J)*DERIV(K)        
 
405
  130         CONTINUE                                                         
 
406
  140     CONTINUE                                                             
 
407
  150 CONTINUE                                                                 
 
408
      DO 170 J = 1,NTERMS                                                      
 
409
          DO 160 K = 1,J                                                       
 
410
              ALPHA(K,J) = ALPHA(J,K)                                          
 
411
  160     CONTINUE                                                             
 
412
  170 CONTINUE                                                                 
 
413
C                                                                              
 
414
C evaluate chi square at starting point                                        
 
415
C                                                                              
 
416
      DO 180 I = 1,NPTS                                                        
 
417
          YFIT(I) = FUNCT(X,I,A)                                               
 
418
  180 CONTINUE                                                                 
 
419
      CHISQ1 = CHISF(Y,SIGMAY,NPTS,NFREE,MODE,YFIT)                            
 
420
C                                                                              
 
421
C invert matrix                                                                
 
422
C                                                                              
 
423
  190 DO 210 J = 1,NTERMS                                                      
 
424
          DO 200 K = 1,NTERMS                                                  
 
425
              IF (ABS(ALPHA(J,J)).LT.1.E-30 .OR.
 
426
     +                ABS(ALPHA(K,K)).LT.1.E-30) THEN
 
427
                 CALL STTPUT
 
428
     +                ('*** WARNING: Insufficient accuracy: NO RESULT',
 
429
     +                ISTAT)
 
430
                 CALL STTPUT
 
431
     +                ('              Scale your input data first', 
 
432
     +                ISTAT)
 
433
                 GOTO 987
 
434
              ENDIF
 
435
              ARRAY(J,K) = ALPHA(J,K)/SQRT(ALPHA(J,J)*ALPHA(K,K))              
 
436
  200     CONTINUE                                                             
 
437
          ARRAY(J,J) = 1. + FLAMDA                                             
 
438
  210 CONTINUE                                                                 
 
439
      CALL INVMAT(ARRAY,NTERMS,DET)                                            
 
440
      DO 230 J = 1,NTERMS                                                      
 
441
          B(J)   = A(J)                                                        
 
442
          DO 220 K = 1,NTERMS                                                  
 
443
              B(J)   = B(J) + BETA(K)*ARRAY(J,K)/                              
 
444
     +                 SQRT(ALPHA(J,J)*ALPHA(K,K))                             
 
445
  220     CONTINUE                                                             
 
446
  230 CONTINUE                                                                 
 
447
C                                                                              
 
448
C if chi square increased, increase FLAMDA and try again                       
 
449
C                                                                              
 
450
      DO 240 I = 1,NPTS                                                        
 
451
          YFIT(I) = FUNCT(X,I,B)                                               
 
452
  240 CONTINUE                                                                 
 
453
      CHISQR = CHISF(Y,SIGMAY,NPTS,NFREE,MODE,YFIT)                            
 
454
      IF (CHISQ1-CHISQR) 250,270,270                                           
 
455
C                                                                              
 
456
C check loops                                                                  
 
457
C                                                                              
 
458
  250 ICNT   = ICNT + 1                                                        
 
459
      IF (ICNT.LT.60) GO TO 260                                                
 
460
  987 IERR   = 1                                                               
 
461
      RETURN                                                                   
 
462
                                                                               
 
463
  260 FLAMDA = 10.*FLAMDA                                                      
 
464
      GO TO 190                                                                
 
465
C                                                                              
 
466
C evaluate parameters                                                          
 
467
C                                                                              
 
468
  270 DO 280 J = 1,NTERMS                                                      
 
469
          A(J)   = B(J)                                                        
 
470
  280 CONTINUE                                                                 
 
471
      FLAMDA = FLAMDA/10.                                                      
 
472
      RETURN                                                                   
 
473
                                                                               
 
474
      END                                                                      
 
475
 
 
476
      DOUBLE PRECISION FUNCTION FUNCT(X,I,A)
 
477
C
 
478
C evaluate terms of function for non-linear least-squares search
 
479
C with form of gaussian peak
 
480
C
 
481
C PARAMETERS:
 
482
C X array of data points
 
483
C I index of data points
 
484
C A array of parameters
 
485
C
 
486
      INTEGER I
 
487
      DOUBLE PRECISION X(*),A(*),Z,XI,Z2
 
488
 
 
489
      XI     = X(I)
 
490
      FUNCT  = A(4)
 
491
      Z      = (XI-A(2))/A(3)
 
492
      Z2     = Z*Z
 
493
      IF (Z2-50.) 40,50,50
 
494
   40 FUNCT  = FUNCT + A(1)*EXP(-Z2*0.5)
 
495
   50 RETURN
 
496
 
 
497
      END
 
498
 
 
499
      DOUBLE PRECISION FUNCTION CHISF(Y,SIGMAY,NPTS,NFREE,MODE,YFIT)
 
500
C
 
501
C evaluate reduced chi square for fit to data
 
502
C chisf = sum((y-yfit)**2/sigma**2)/nfree
 
503
C
 
504
C PARAMETERS:
 
505
C Y array of data points
 
506
C SIGMAY array of standard deviations
 
507
C NPTS number of data points
 
508
C NFREE number of degrees of freedom
 
509
C MODE determines method of weighting least-squares fit
 
510
C  +1 (instrumental) weight = 1./SIGMAY(i)**2
 
511
C   0 (no weighting) weight = 1.
 
512
C  -1 (statistical ) weight = 1./Y(i)
 
513
C YFIT array of calculated values of Y
 
514
C
 
515
      INTEGER NFREE,MODE,NPTS,I
 
516
      REAL             Y(1)
 
517
CC      DOUBLE PRECISION SIGMAY(1),YFIT(1),CHISF
 
518
      DOUBLE PRECISION SIGMAY(1),YFIT(1)
 
519
      DOUBLE PRECISION CHISQ,WEIGHT
 
520
C
 
521
      CHISQ  = 0.
 
522
      IF (NFREE) 30,30,40
 
523
   30 CHISF  = 0.
 
524
      RETURN
 
525
C
 
526
C accumulate chi square
 
527
C
 
528
   40 DO 110 I = 1,NPTS
 
529
          IF (MODE) 50,80,90
 
530
   50     IF (Y(I)) 70,80,60
 
531
   60     WEIGHT = 1./Y(I)
 
532
          GO TO 100
 
533
 
 
534
   70     WEIGHT = 1./ (-Y(I))
 
535
          GO TO 100
 
536
 
 
537
   80     WEIGHT = 1.
 
538
          GO TO 100
 
539
 
 
540
   90     WEIGHT = 1./SIGMAY(I)**2
 
541
  100     CHISQ  = CHISQ + WEIGHT* (Y(I)-YFIT(I))**2
 
542
  110 CONTINUE
 
543
C
 
544
C divide by number of degrees of freedom
 
545
C
 
546
      CHISF  = CHISQ/NFREE
 
547
      RETURN
 
548
 
 
549
      END
 
550
 
 
551
      SUBROUTINE FDERI(X,I,A,DELTAA,NTERMS,DERIV)                              
 
552
C                                                                              
 
553
C evaluates derivatives of function for least squares search                   
 
554
C with form of gaussian peak                                                   
 
555
C                                                                              
 
556
C PARAMETERS:                                                                  
 
557
C X array of data points of independent variable                               
 
558
C I index of data points                                                       
 
559
C A array of parameters                                                        
 
560
C DELTAA array of parameter increments                                         
 
561
C NTERMS number of parameters                                                  
 
562
C DERIV derivatives of function                                                
 
563
C                                                                              
 
564
      INTEGER I,J,NTERMS
 
565
 
 
566
      DOUBLE PRECISION X(*),A(*),DELTAA(*),DERIV(*) 
 
567
      DOUBLE PRECISION Z,Z2,XI
 
568
                                                                               
 
569
      XI     = X(I)                                                            
 
570
      Z      = (XI-A(2))/A(3)                                                  
 
571
      Z2     = Z*Z                                                             
 
572
      IF (Z2-50.) 40,20,20                                                     
 
573
   20 DO 30 J = 1,3                                                            
 
574
          DERIV(J) = 0.                                                        
 
575
   30 CONTINUE                                                                 
 
576
      GO TO 50                                                                 
 
577
C                                                                              
 
578
C analytical expression for derivatives                                        
 
579
C                                                                              
 
580
   40 DERIV(1) = EXP(-Z2/2.)                                                   
 
581
      DERIV(2) = A(1)*DERIV(1)*Z/A(3)                                          
 
582
      DERIV(3) = DERIV(2)*Z                                                    
 
583
   50 DERIV(4) = 1.                                                            
 
584
      RETURN                                                                   
 
585
                                                                               
 
586
      END                                                                      
 
587
 
 
588
      SUBROUTINE INVMAT(ARR,N,DET)
 
589
C                                         all rights reserved
 
590
C
 
591
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  0:26 - 4 DEC 1987
 
592
C
 
593
C.LANGUAGE: F77+ESOext
 
594
C
 
595
C.AUTHOR: D.PONZ
 
596
C         inverts a symmetric matrix and calculates its determinant
 
597
C         using a gauss-jordan elimination method. the degree n of
 
598
C         the matrix must not be greater than 10.
 
599
C         the inversion is in double precision, and
 
600
C         the array ARR must be double. the determinant is allways real.
 
601
C
 
602
C   PARAMS :
 
603
C         ARR :  array containing the (n*n) matrix to be inverted.
 
604
C                the inverted matrix is returned in ARR.
 
605
C         N   :  degree of matrix
 
606
C         DET :  determinant of input-matrix
 
607
C
 
608
C      IMPLICIT NONE
 
609
                                         
 
610
      DOUBLE PRECISION DET
 
611
      DOUBLE PRECISION ARR(10,10),AMAX,SAVE,DDET
 
612
      INTEGER IK(10),JK(10),I,J,K,N,L
 
613
C
 
614
      DET    = 1.0
 
615
      DDET   = 1.0D0
 
616
      DO 210 K = 1,N
 
617
          AMAX   = 0.0D0
 
618
C
 
619
C  find largest element ARR(i,j) in rest of matrix
 
620
C
 
621
   10     DO 40 I = K,N
 
622
              DO 30 J = K,N
 
623
                  IF (DABS(AMAX)-DABS(ARR(I,J))) 20,20,30
 
624
   20             AMAX   = ARR(I,J)
 
625
                  IK(K)  = I
 
626
                  JK(K)  = J
 
627
   30         CONTINUE
 
628
   40     CONTINUE
 
629
C  ----------------------------------------------------
 
630
C  interchange rows and columns to put AMAX in ARR(k,k)
 
631
C  ----------------------------------------------------
 
632
          IF (AMAX) 60,50,60
 
633
   50     DET    = 0.
 
634
          GO TO 280
 
635
 
 
636
   60     I      = IK(K)
 
637
          IF (I-K) 10,90,70
 
638
   70     DO 80 J = 1,N
 
639
              SAVE   = ARR(K,J)
 
640
              ARR(K,J) = ARR(I,J)
 
641
              ARR(I,J) = -SAVE
 
642
   80     CONTINUE
 
643
C
 
644
   90     J      = JK(K)
 
645
          IF (J-K) 10,120,100
 
646
  100     DO 110 I = 1,N
 
647
              SAVE   = ARR(I,K)
 
648
              ARR(I,K) = ARR(I,J)
 
649
              ARR(I,J) = -SAVE
 
650
  110     CONTINUE
 
651
C  -------------------------------------
 
652
C  accumulate elements of inverse matrix
 
653
C  -------------------------------------
 
654
  120     DO 140 I = 1,N
 
655
              IF (I-K) 130,140,130
 
656
  130         ARR(I,K) = -ARR(I,K)/AMAX
 
657
  140     CONTINUE
 
658
C
 
659
          DO 180 I = 1,N
 
660
              DO 170 J = 1,N
 
661
                  IF (I-K) 150,170,150
 
662
  150             IF (J-K) 160,170,160
 
663
  160             ARR(I,J) = ARR(I,J) + ARR(I,K)*ARR(K,J)
 
664
  170         CONTINUE
 
665
  180     CONTINUE
 
666
C
 
667
          DO 200 J = 1,N
 
668
              IF (J-K) 190,200,190
 
669
  190         ARR(K,J) = ARR(K,J)/AMAX
 
670
  200     CONTINUE
 
671
          ARR(K,K) = 1.0/AMAX
 
672
          DET    = DET*AMAX
 
673
          DDET   = DDET*AMAX
 
674
C
 
675
  210 CONTINUE
 
676
C  --------------------------
 
677
C  restore ordering of matrix
 
678
C  --------------------------
 
679
      DO 270 L = 1,N
 
680
          K      = N - L + 1
 
681
          J      = IK(K)
 
682
          IF (J-K) 240,240,220
 
683
  220     DO 230 I = 1,N
 
684
              SAVE   = ARR(I,K)
 
685
              ARR(I,K) = -ARR(I,J)
 
686
              ARR(I,J) = SAVE
 
687
  230     CONTINUE
 
688
C
 
689
  240     I      = JK(K)
 
690
          IF (I-K) 270,270,250
 
691
  250     DO 260 J = 1,N
 
692
              SAVE   = ARR(K,J)
 
693
              ARR(K,J) = -ARR(I,J)
 
694
              ARR(I,J) = SAVE
 
695
  260     CONTINUE
 
696
C
 
697
  270 CONTINUE
 
698
      DET    = DDET
 
699
      RETURN
 
700
C  -------------------
 
701
C  determinant is zero
 
702
C  -------------------
 
703
  280 DET    = 0.0
 
704
      RETURN
 
705
 
 
706
      END