← Back to branch summary

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

« back to all changes in this revision

Viewing changes to stdred/echelle/src/necripp1.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:
stdred/echelle/src/necripp1.for
 
1
C @(#)necripp1.for      19.1 (ESO-DMD) 02/25/03 14:20:25
 
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 Massachusetts Ave, Cambridge, 
 
18
C MA 02139, USA.
 
19
C
 
20
C Correspondence 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
       PROGRAM ECHRP1
 
30
C
 
31
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
32
C
 
33
C.IDENTIFICATION
 
34
C
 
35
C  program  ECHRIPP1              version 1.0       870820
 
36
C
 
37
C  Author :   J.D.Ponz            ESO - Garching
 
38
C
 
39
C.MODIFICATIONS
 
40
C
 
41
C  version 1.1   P. Ballester   901009   Option to freeze K and Alpha
 
42
C
 
43
C.KEYWORDS
 
44
C
 
45
C  ECHELLE, CASPEC, BLAZE FUNCTION
 
46
C
 
47
C.PURPOSE
 
48
C
 
49
C  COMPUTE THE ECHELLE CONSTANTS ALPHA AND K FOR EACH ECHELLE ORDER
 
50
C  (only if the flag is not set to FREEZE)
 
51
C
 
52
C.ALGORITHM
 
53
C
 
54
C  a new one...
 
55
C
 
56
C.INPUT/OUTPUT
 
57
C
 
58
C  RIPPLE/ECHELLE INPUT OUTPUT  K0,A0
 
59
C
 
60
C.VERSION
 
61
C 
 
62
C 010713                last modif
 
63
C 
 
64
C-------------------------------------------------------------
 
65
C 
 
66
       IMPLICIT  NONE
 
67
C 
 
68
       INTEGER          NPIXA(2),NPTOT(100),NORDER(100)
 
69
       INTEGER          I, ISTAT, KUN, KNUL, INIMA, OUTIMA
 
70
       INTEGER          MADRID(1), NAXISA
 
71
       INTEGER*8        IPNTRA, IPNTRB
 
72
C 
 
73
       DOUBLE PRECISION WSTART(100),STEPA(2),STARTA(2)
 
74
       DOUBLE PRECISION K0, ALPHA0,LAMBDA1,LAMBDA2
 
75
C 
 
76
       REAL             CUT(4), PARS(4), XMIN, XMAX
 
77
C 
 
78
       CHARACTER*80     INFRM, OUTFRM, TABLE, FREEZE
 
79
       CHARACTER        IDENTA*72,CUNITA*64
 
80
C 
 
81
       INCLUDE          'MID_INCLUDE:ST_DEF.INC'
 
82
       COMMON/VMR/MADRID
 
83
       INCLUDE          'MID_INCLUDE:ST_DAT.INC'
 
84
C
 
85
C ... INITIALIZE SYSTEM
 
86
C
 
87
       CALL STSPRO('ECHRP1')
 
88
       CALL STKRDC('P1',1,1,80,I,INFRM,KUN,KNUL,ISTAT)
 
89
       CALL STKRDC('P2',1,1,80,I,OUTFRM,KUN,KNUL,ISTAT)
 
90
       CALL CLNFRA(INFRM,INFRM,0)
 
91
       CALL CLNFRA(OUTFRM,OUTFRM,0)
 
92
       CALL STKRDR('INPUTR',1,4,I,PARS,KUN,KNUL,ISTAT)
 
93
       TABLE = 'riptab'
 
94
C       CALL STKRDC('P4',1,1,80,I,TABLE,KUN,KNUL,ISTAT)
 
95
       CALL STKRDC ('P5',1,1,80,I,FREEZE,KUN,KNUL,ISTAT)
 
96
       CALL STTPUT(' Ripple correction. Method FIT',ISTAT)
 
97
       LAMBDA1 = PARS(1)
 
98
       LAMBDA2 = PARS(2)
 
99
       K0     = PARS(3)
 
100
       ALPHA0 = PARS(4)
 
101
C
 
102
C ... MAP INPUT/OUTPUT FRMS
 
103
C
 
104
       CALL STIGET(INFRM,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
 
105
     .             2,NAXISA,NPIXA,STARTA,STEPA,IDENTA,
 
106
     .             CUNITA,IPNTRA,INIMA,ISTAT)
 
107
       CALL STDRDR(INIMA,'LHCUTS',1,4,I,CUT,KUN,KNUL,ISTAT)
 
108
       CALL STDRDD(INIMA,'WSTART',1,NPIXA(2),I,WSTART,KUN,KNUL,ISTAT)
 
109
       CALL STDRDI(INIMA,'NPTOT',1,NPIXA(2),I,NPTOT,KUN,KNUL,ISTAT)
 
110
       CALL STDRDI(INIMA,'NORDER',1,NPIXA(2),I,NORDER,KUN,KNUL,ISTAT)
 
111
       CALL STIPUT(OUTFRM,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
 
112
     .             NAXISA,NPIXA,STARTA,STEPA,IDENTA,
 
113
     .             CUNITA,IPNTRB,OUTIMA,ISTAT)
 
114
C
 
115
C ... RIPPLE FUNCTION
 
116
C
 
117
       CALL RIPPLE1(NPIXA(1),NPIXA(2),MADRID(IPNTRA),MADRID(IPNTRB),
 
118
     .      WSTART,STEPA(1),NPTOT,NORDER,K0,ALPHA0,XMIN,XMAX,TABLE,
 
119
     .      LAMBDA1,LAMBDA2,FREEZE)
 
120
       XMIN = CUT(3)
 
121
       XMIN = CUT(4)
 
122
C
 
123
C ... WRITE PROCESS DESCRIPTORS
 
124
C
 
125
       CUT(1) = XMIN
 
126
       CUT(2) = XMAX
 
127
       CUT(3) = XMIN
 
128
       CUT(4) = XMAX
 
129
       CALL DSCUPT(INIMA,OUTIMA,' ',ISTAT)
 
130
       CALL STDWRR(OUTIMA,'LHCUTS',CUT,1,4,KUN,ISTAT)
 
131
       CALL STSEPI
 
132
       END
 
133
 
 
134
       SUBROUTINE RIPPLE1(NX,NY,X,Y,WS,WSTEP,NP,NO,K0,A0,
 
135
     .                        XMIN,XMAX,TABLE,LAMBDA1,LAMBDA2,FREEZE)
 
136
C
 
137
C       CORRECT FOR BLAZE EFFECT USING THE OVERLAPPING REGION
 
138
C
 
139
C       INPUT ARGUMENTS
 
140
C       NX       INTG       MAXIMUM NUMBER OF SAMPLES PER ORDER
 
141
C       NY       INTG       NUMBER OF ORDERS
 
142
C       X       REAL       X(NX,NY) INPUT IMAGE
 
143
C       WS       REAL       STARTING WAVELENGTH PER EACH ORDER
 
144
C       WSTEP       REAL       WAVELENGTH STEP
 
145
C       NP       INTG       NUMBER OF USEFUL SAMPLES PER ORDER
 
146
C       NO       INTG       SPECTRAL ORDER NUMBER
 
147
C       K0       REAL       INITIAL GUESS FOR THE CONSTANT K
 
148
C       A0       REAL       INITIAL GUESS FOR ALPHA
 
149
C       FREEZE   CHARACTER  Flag (NULL/FREEZE) to freeze values of K and Alpha
 
150
C
 
151
C       OUTPUT ARGUMENTS
 
152
C       Y       REAL       Y(NX,NY) CORRECTED IMAGE
 
153
C       XMIN       REAL       MINIMUN AND
 
154
C       XMAX       REAL       MAXIMUM VALUES
 
155
C
 
156
       IMPLICIT NONE
 
157
C 
 
158
       INTEGER           NX, NY
 
159
       INTEGER           NP(NY),NO(NY),IERR(2),ICOL(7), TID
 
160
       INTEGER           ISTAT, ISEQ, I, IORDER, I1, I2, ILR
 
161
C 
 
162
       REAL              X(NX,NY), Y(NX,NY), XMIN, XMAX
 
163
       REAL              VAL(7)
 
164
C 
 
165
       DOUBLE PRECISION  WS(NY), WSTEP, K0, A0,LAMBDA1,LAMBDA2
 
166
       DOUBLE PRECISION  K(2), KM, AM, A(2)
 
167
       DOUBLE PRECISION  DVAL1, DVAL2
 
168
C 
 
169
       LOGICAL           NULL(2)
 
170
C 
 
171
       CHARACTER*16      LABEL(7),UNIT(7),FORM
 
172
       CHARACTER*80      LINE,FREEZE
 
173
       CHARACTER*(*)     TABLE
 
174
C 
 
175
       INCLUDE          'MID_INCLUDE:ST_DEF.INC'
 
176
       INCLUDE          'MID_INCLUDE:ST_DAT.INC'
 
177
       DATA              FORM/'G12.4'/
 
178
       DATA              LABEL(1)/'ORDER'/,UNIT(1)/'UNITLESS'/
 
179
       DATA              LABEL(2)/'K'/,UNIT(2)/'UNITLESS'/
 
180
       DATA              LABEL(3)/'ALPHA'/,UNIT(3)/'UNITLESS'/
 
181
       DATA              LABEL(4)/'KFIT'/,UNIT(4)/'UNITLESS'/
 
182
       DATA              LABEL(5)/'AFIT'/,UNIT(5)/'UNITLESS'/
 
183
       DATA              LABEL(6)/'YMAX'/,UNIT(6)/'FLUX'/
 
184
       DATA              LABEL(7)/'STATUS'/,UNIT(7)/'UNITLESS'/
 
185
       DATA              TID /-1/
 
186
C
 
187
       ISEQ = 0
 
188
       XMIN = 1.E30
 
189
       XMAX = -XMIN
 
190
       IF (K0.LE.0.0) THEN                     
 
191
C
 
192
C ... CORRECT RIPPLE IF K0 IS NEGATIVE
 
193
C
 
194
         CALL TBTOPN(TABLE,F_I_MODE,TID,ISTAT)
 
195
         CALL TBLSER(TID,LABEL(2),ICOL(1),ISTAT)
 
196
         CALL TBLSER(TID,LABEL(3),ICOL(2),ISTAT)
 
197
         DO 10 IORDER=1, NY
 
198
           CALL TBRRDR(TID,IORDER,2,ICOL,VAL,NULL,ISTAT)
 
199
           DVAL1 = VAL(1)
 
200
           DVAL2 = VAL(2)
 
201
           CALL RIPCOR(WS(IORDER),WSTEP,NP(IORDER),X(1,IORDER),
 
202
     .                 Y(1,IORDER),NO(IORDER),DVAL1,DVAL2,NX)
 
203
10       CONTINUE
 
204
       ELSE
 
205
C
 
206
C ... FIND BLAZE PARAMETERS
 
207
C
 
208
         CALL TBTINI(TABLE,F_TRANS,F_O_MODE,8,NY,TID,ISTAT)
 
209
         DO 20 I = 1, 6
 
210
           CALL TBCINI(TID,D_R4_FORMAT,1,FORM,UNIT(I),LABEL(I),
 
211
     .                 ICOL(I),ISTAT)
 
212
20       CONTINUE
 
213
C
 
214
C ... ITERATION ON ALL CONSECUTIVE PAIRS OF ORDERS
 
215
C
 
216
         CALL STTPUT(' ORDER NUMER  FITTED K  FITTED ALPHA',ISTAT)
 
217
         CALL STTPUT(' ----------- ---------- ------------',ISTAT)
 
218
         DO 40 IORDER = 2, NY-1
 
219
 
 
220
         IF (FREEZE(1:1).EQ.'F'.OR.FREEZE(1:1).EQ.'f') THEN
 
221
 
 
222
           KM = K0
 
223
           AM = A0
 
224
 
 
225
         ELSE
 
226
 
 
227
           I1 = IORDER-2
 
228
C
 
229
C ... CONSIDER LEFT AND RIGHT OVERLAPPING REGIONS
 
230
C
 
231
           DO 30 ILR = 0, 1
 
232
             I1 = I1 + 1
 
233
             I2 = I1 + 1
 
234
             CALL FITKA(WS(I1),WSTEP,NP(I1),X(1,I1),NO(I1),     
 
235
     .                  WS(I2),NP(I2),X(1,I2),NO(I2),
 
236
     .                  K0,A0,K(ILR+1),A(ILR+1),IERR(ILR+1),
 
237
     .                  LAMBDA1,LAMBDA2)
 
238
30         CONTINUE
 
239
C
 
240
C ... COMPUTE MEAN VALUE
 
241
C
 
242
           KM = 0.5*(K(1) + K(2))
 
243
           AM = 0.5*(A(1) + A(2))
 
244
 
 
245
         ENDIF          ! End of test on FREEZE
 
246
 
 
247
C           TYPE *,IORDER,KM,AM
 
248
C
 
249
C ... CORRECT FOR RIPPLE
 
250
C
 
251
           CALL RIPCOR(WS(IORDER),WSTEP,NP(IORDER),X(1,IORDER),
 
252
     .                     Y(1,IORDER),NO(IORDER),KM,AM,NX)
 
253
           IF (IORDER.EQ.2) THEN
 
254
             CALL RIPCOR(WS(1),WSTEP,NP(1),X(1,1),Y(1,1),
 
255
     .                       NO(1),KM,AM,NX)
 
256
             VAL(1) = NO(1)
 
257
             VAL(2) = KM
 
258
             VAL(3) = AM
 
259
             VAL(4) = KM
 
260
             VAL(5) = AM
 
261
             ISEQ = ISEQ + 1
 
262
             CALL TBRWRR(TID,ISEQ,6,ICOL,VAL,ISTAT)
 
263
           ENDIF
 
264
           VAL(1) = NO(IORDER)
 
265
           VAL(2) = KM
 
266
           VAL(3) = AM
 
267
           VAL(4) = KM
 
268
           VAL(5) = AM
 
269
           ISEQ = ISEQ + 1
 
270
           CALL TBRWRR(TID,ISEQ,6,ICOL,VAL,ISTAT)
 
271
           IF (IORDER.EQ.NY-1) THEN
 
272
             CALL RIPCOR(WS(NY),WSTEP,NP(NY),X(1,NY),Y(1,NY),
 
273
     .                   NO(NY),KM,AM,NX)
 
274
             VAL(1) = NO(NY)
 
275
             VAL(2) = KM
 
276
             VAL(3) = AM
 
277
             VAL(4) = KM
 
278
             VAL(5) = AM
 
279
             ISEQ = ISEQ + 1
 
280
             CALL TBRWRR(TID,ISEQ,6,ICOL,VAL,ISTAT)
 
281
           ENDIF
 
282
C
 
283
C ... WRITE OUT RESULTS
 
284
C
 
285
           WRITE(LINE,100) NO(IORDER),KM,AM
 
286
           CALL STTPUT(LINE,ISTAT)
 
287
40       CONTINUE
 
288
         CALL TBSINI(TID,ISTAT)
 
289
       ENDIF
 
290
       CALL TBTCLO(TID,ISTAT)
 
291
       RETURN
 
292
100    FORMAT(4X,I8,F12.2,F12.4)
 
293
       END
 
294
 
 
295
       SUBROUTINE RIPCOR(W1,WS,NP,X,Y,M,K,A,NX)
 
296
C
 
297
C       CORRECT FOR RIPPLE FUNCTION
 
298
C
 
299
C       INPUT ARGUMENTS
 
300
C       W1      INITIAL WAVELENGTH
 
301
C       WS      WAVELENGTH STEP
 
302
C       NP      NO. OF POINTS
 
303
C       X       INPUT FLUX
 
304
C       M       ORDER NUMBER
 
305
C       K       K CONSTANT
 
306
C       A       ALPHA CONSTANT
 
307
C
 
308
C       OUTPUT ARGUMENTS
 
309
C       Y       CORRECTED FLUX
 
310
C
 
311
       IMPLICIT NONE
 
312
       INTEGER          NP, NX, M
 
313
       REAL             X(NP),Y(NX)
 
314
       DOUBLE PRECISION W1,WS,K,A
 
315
       DOUBLE PRECISION PI,DK,DA,DW,DM,R,DX,PA,DC
 
316
       INTEGER          I
 
317
       DATA             PI/3.1415926535897932D0/
 
318
C
 
319
       DK = K
 
320
       DA = A
 
321
       DM = M
 
322
       PA = PI*A
 
323
       DC = DM/DK
 
324
       DO 10 I = 1, NP
 
325
         DW   = W1 + (I-1)*WS
 
326
         DX   = (PA*DM*DC)*(DW-1.D0/DC)
 
327
         IF (DABS(DX).LT.1.D-10) THEN
 
328
           Y(I) = X(I)
 
329
         ELSE                  
 
330
           R    = DSIN(DX)/DX
 
331
           Y(I) = X(I)/(R*R) 
 
332
         ENDIF
 
333
10     CONTINUE
 
334
       DO 20 I = NP+1, NX
 
335
         Y(I) = 0.
 
336
20     CONTINUE
 
337
       RETURN
 
338
       END
 
339
 
 
340
       SUBROUTINE FITKA(WSTART1,WSTEP,N1,X1,NO1,     
 
341
     .                  WSTART2,N2,X2,NO2,K0,A0,K,A,IERR,
 
342
     .                  LAMBDA1,LAMBDA2)
 
343
C
 
344
C       COMPUTE CONTANT K AND ALPHA
 
345
C
 
346
C       INPUT ARGUMENTS
 
347
C       WSTEP    WAVELENGTH STEP
 
348
C       WSTART1  STARTING WAVELENGTH ORDER M+1
 
349
C       N1       NUMBER OF POINTS IN ORDER M+1
 
350
C       X1       FLUX IN ORDER M+1
 
351
C       NO1      ORDER M+1
 
352
C       WSTART2  STARTING WAVELENGTH ORDER M
 
353
C       N2       NUMBER OF POINTS IN ORDER M
 
354
C       X2       FLUX IN ORDER M
 
355
C       NO2      ORDER M
 
356
C       K0       INITIAL K VALUE
 
357
C       A0       ALPHA
 
358
C
 
359
C       OUTPUT ARGUMENTS
 
360
C       K       FINAL K VALUE
 
361
C       A       FINAL A VALUE
 
362
C       IERR    CONVERGENCE FLAG 0 - GOOD VALUE
 
363
C                               OTHERWISE IT DOES NOT CONVERGE
 
364
C
 
365
C
 
366
       IMPLICIT NONE
 
367
       INTEGER           N1, N2
 
368
       REAL              X1(N1), X2(N2)
 
369
       DOUBLE PRECISION  WSTART1, WSTART2, WSTEP
 
370
       DOUBLE PRECISION  DW1, DW2, W2, K0, K, A0, A,LAMBDA1,LAMBDA2
 
371
       DOUBLE PRECISION  ETA, XTOL, STEPMX, XPAR(2), FSUMSQ
 
372
       DOUBLE PRECISION  FVEC(300), FJAC(300,2),S(2), V(2,2), W(1500)
 
373
       DOUBLE PRECISION  XX1(300),XX2(300),XSTART,XSTEP
 
374
       INTEGER           ORD1, ORD2, IW(1), ISTAT, NPIX, I1, I2
 
375
       INTEGER           NO1, NO2, J, J1, IPRINT, NPAR, MAXCAL
 
376
       INTEGER           IERR, NF, NITER, IFAIL, LV, LW, LJ, LIW
 
377
C       INTEGER           I
 
378
C       CHARACTER*80      LINE
 
379
C
 
380
       DOUBLE PRECISION  DSQRT, X02AAF
 
381
       EXTERNAL          LSQFUN, LSQMON
 
382
       COMMON/E04PAR/XSTART,XSTEP,ORD1,ORD2,XX1,XX2
 
383
C
 
384
       DW1  = (WSTART1 + (N1-1)*WSTEP) - WSTART2
 
385
       IF (DW1.LE.0.D0) THEN
 
386
         CALL STTPUT('Warning: There is no order overlap',ISTAT)
 
387
         RETURN
 
388
       ENDIF
 
389
C       DW2  = DW1*0.8D0
 
390
C       W2   = DW1*0.1D0
 
391
       DW2 = LAMBDA2
 
392
       W2 = LAMBDA1
 
393
       IF ((DW2+W2).GT.DW1) THEN 
 
394
          CALL STTPUT('Warning: Wrong wavelengths !',ISTAT)
 
395
C          TYPE *,DW1
 
396
       ENDIF
 
397
C       TYPE *,' *** ORDERS :',NO1,NO2,' ***',W2,DW2
 
398
       NPIX = DW2/WSTEP
 
399
       NPIX = MIN(NPIX,300)
 
400
       I1   = (WSTART2+W2-WSTART1)/WSTEP 
 
401
       I2   = W2/WSTEP 
 
402
C
 
403
C ... PREPARE INPUT
 
404
C
 
405
       ORD1 = NO1
 
406
       ORD2 = NO2
 
407
       XSTEP  = WSTEP
 
408
       XSTART = WSTART2 + W2 
 
409
       DO 10 J = 1, NPIX
 
410
         J1 = J-1
 
411
         XX1(J) = X1(I1 + J1)
 
412
         XX2(J) = X2(I2 + J1)
 
413
10     CONTINUE
 
414
C100    CONTINUE
 
415
C
 
416
C ... START MINIMIZATION 
 
417
C
 
418
       IPRINT = 1
 
419
       NPAR   = 2
 
420
       MAXCAL = 50*NPAR
 
421
       ETA    = 0.9
 
422
       XTOL   = 10.*DSQRT(X02AAF(XTOL))
 
423
       STEPMX = 10.
 
424
       XPAR(1)= K0
 
425
       XPAR(2)= A0
 
426
       LV     = 2
 
427
       LJ     = 300
 
428
       LIW    = 1
 
429
       LW     = 1500
 
430
       IFAIL  = 1
 
431
       CALL E04GDF(NPIX,NPAR,LSQFUN,LSQMON,IPRINT,MAXCAL,ETA,XTOL,
 
432
     .      STEPMX,XPAR,FSUMSQ,FVEC,FJAC,LJ,S,V,LV,NITER,NF,
 
433
     .      IW,LIW,W,LW,IFAIL)
 
434
C
 
435
C ... TYPE RESULTS
 
436
C
 
437
C       TYPE *,'IFAIL,SUMSQ',IFAIL,FSUMSQ
 
438
C       TYPE *,'K,A',XPAR(1),XPAR(2)
 
439
       K = XPAR(1)
 
440
       A = XPAR(2)
 
441
       RETURN
 
442
       END
 
443
 
 
444
       SUBROUTINE LSQFUN(IFLAG,M,N,XC,FVECC,FJACC,LJC,IW,LIW,W,LW)
 
445
C
 
446
C       SUBROUTINE TO COMPUTE THE RESIDUALS AND THEIR FIRST DERIVATIVE
 
447
C
 
448
C
 
449
       IMPLICIT NONE
 
450
       INTEGER              IFLAG, M, N, LJC, LIW, LW
 
451
C
 
452
       DOUBLE PRECISION     FJACC(LJC,N), FVECC(M), W(LW), XC(N)       
 
453
C
 
454
       DOUBLE PRECISION     XX1(300),XX2(300),XSTART,XSTEP
 
455
       INTEGER              ORD1, ORD2, IW(LIW), I
 
456
C
 
457
       DOUBLE PRECISION     PA, PI, WL, PA2, P2
 
458
       DOUBLE PRECISION     X10, X1, X13, S1, C1, R1, DKR1, DAR1
 
459
       DOUBLE PRECISION     X20, X2, X23, S2, C2, R2, DKR2, DAR2
 
460
C
 
461
       COMMON/E04PAR/XSTART,XSTEP,ORD1,ORD2,XX1,XX2
 
462
       DATA              PI/3.1415926535897932D0/
 
463
C
 
464
       PA = PI*XC(2)
 
465
       PA2 = 2.D0*PA
 
466
       P2  = 2.D0*PI
 
467
       DO 10 I = 1, M
 
468
         WL  = XSTART + (I-1)*XSTEP
 
469
         X10 = ORD1 - XC(1)/WL
 
470
         X1  = PA*X10
 
471
         X13 = X1*X1*X1
 
472
         S1  = DSIN(X1)
 
473
         C1  = DCOS(X1)
 
474
         R1  = S1/X1
 
475
         X20 = ORD2 - XC(1)/WL
 
476
         X2  = PA*X20
 
477
         X23 = X2*X2*X2
 
478
         S2  = DSIN(X2)
 
479
         C2  = DCOS(X2)
 
480
         R2  = S2/X2
 
481
         IF (IFLAG.EQ.2) FVECC(I) = R1*R1/XX1(I) - R2*R2/XX2(I)
 
482
         DKR1 = (PA2/(WL*X13))*(S1*S1-X1*S1*C1) 
 
483
         DKR2 = (PA2/(WL*X23))*(S2*S2-X2*S2*C2)
 
484
         DAR1 = (X1*S1*C1-S1*S1)*P2*X10/X13
 
485
         DAR2 = (X2*S2*C2-S2*S2)*P2*X20/X23
 
486
         FJACC(I,1) = DKR1/XX1(I) - DKR2/XX2(I)
 
487
         FJACC(I,2) = DAR1/XX1(I) - DAR2/XX2(I)
 
488
10     CONTINUE
 
489
       RETURN
 
490
       END
 
491
 
 
492
       SUBROUTINE LSQMON(M,N,XC,FVECC,FJACC,LJC,S,IGRADE,NITER,NF,
 
493
     .                       IW,LIW,W,LW)
 
494
C
 
495
C       MONITORING SUBROUTINE
 
496
C
 
497
       IMPLICIT NONE
 
498
       INTEGER          M, N, LJC, IGRADE, NITER, NF, LIW, LW
 
499
       DOUBLE PRECISION FJACC(LJC,N), FVECC(M), W(LW), XC(N), S(N)
 
500
       INTEGER          IW(LIW)
 
501
C
 
502
C       DOUBLE PRECISION FSUMSQ, GTG, G(50), F01DEF
 
503
C       INTEGER          J
 
504
C
 
505
C       FSUMSQ = F01DEF(FVECC,FVECC, M)
 
506
C       CALL LSQGRD(M, N, FVECC, FJACC, LJC, G)
 
507
C       GTG    = F01DEF(G,G,N)
 
508
C       WRITE(6, 99999) NITER, NF, FSUMSQ, GTG, IGRADE
 
509
C       WRITE(6, 99998)
 
510
C       DO 20 J = 1,  N
 
511
C          WRITE(6, 99997) XC(J), G(J), S(J)
 
512
C20     CONTINUE
 
513
       RETURN
 
514
C99999  FORMAT(///'    ITNS',4X,'F  EVALS',10X, ' SUMSQ',13X,'GTG',
 
515
C     . 8X,'GRADE'/1X I4,6X,I5,6X,1PE13.5,6X,1PE9.1,6X,I3)
 
516
C99998  FORMAT(/8X,'X',20X,'G',11X,'SINGULAR VALUES')
 
517
C99997  FORMAT(1X,1PE13.5,10X,1PE9.1,10X,1PE9.1)
 
518
       END
 
519
 
 
520
       SUBROUTINE LSQGRD(M, N, FVECC, FJACC, LJC, G)
 
521
C
 
522
       IMPLICIT NONE
 
523
       INTEGER LJC, M, N
 
524
       DOUBLE PRECISION FVECC(M), FJACC(LJC, N), G(N)
 
525
       DOUBLE PRECISION SUM
 
526
       INTEGER          J, I
 
527
C
 
528
       DO 40 J = 1, N
 
529
          SUM = 0.D0
 
530
          DO 20 I =1, M
 
531
             SUM = SUM + FJACC(I, J)*FVECC(I)
 
532
20        CONTINUE
 
533
          G(J) = SUM + SUM
 
534
40     CONTINUE
 
535
C
 
536
       RETURN
 
537
       END