← Back to branch summary

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

« back to all changes in this revision

Viewing changes to contrib/romafot/libsrc/elmor.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/romafot/libsrc/elmor.for
 
1
C @(#)elmor.for 19.1 (ES0-DMD) 02/25/03 13:29:43
 
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
      SUBROUTINE ELMRF(IVX,IVY,VZ,NP,VP,FS,VPES,NC,BETA,SQM,LG,WEI,
 
30
     2                  SIG)
 
31
C+++
 
32
C.PURPOSE: Moffat o Gauss sigma fisso - piano orizzontale
 
33
C---
 
34
      IMPLICIT     NONE
 
35
      INTEGER      NST
 
36
      INTEGER      NIND
 
37
      PARAMETER    (NST=40)
 
38
      PARAMETER    (NIND=4*NST+3)
 
39
C
 
40
      INTEGER      IVX(1), IVY(1)
 
41
      REAL         VZ(1)
 
42
      INTEGER      NP
 
43
      REAL         VP(1)
 
44
      REAL         FS
 
45
      REAL         VPES(7)
 
46
      INTEGER      NC
 
47
      REAL         BETA
 
48
      REAL         SQM
 
49
      INTEGER      LG, THREE
 
50
      REAL         WEI(1)
 
51
      REAL         SIG
 
52
C
 
53
      REAL         VC, VTP, VPI
 
54
      REAL         CO
 
55
      INTEGER      N1
 
56
      INTEGER      INK, INK1, NCD
 
57
      REAL         EPE, EP2, EEE
 
58
      REAL         RMT
 
59
      REAL         VTN, VFZ, VFF, VFU, VG       
 
60
      REAL         DIX, DIY
 
61
      REAL         EP, EPES
 
62
      INTEGER      NIN
 
63
      INTEGER      I, J, K, L, N
 
64
      REAL         FK(NST)
 
65
C
 
66
      COMMON      /SUFR/RMT(NIND,NIND),VTN(NIND),VFZ(NIND),VFF(NIND),
 
67
     2                  VC(NIND)
 
68
C
 
69
C *** start the code
 
70
      THREE  = 3
 
71
      NIN=NC*3
 
72
 
 
73
      DO 10 I=1,NIN
 
74
         VFF(I) = 0
 
75
         VFZ(I) = 0
 
76
         VC(I)  = 0
 
77
         DO 11 J = 1,NIN
 
78
            RMT(I,J) = 0
 
79
   11    CONTINUE
 
80
   10 CONTINUE
 
81
 
 
82
      DO 20 I = 1,NC
 
83
         IF(BETA.LE.0.) THEN
 
84
            FK(I) = -4*ALOG(2.)/(VP(4*I+3)**2)
 
85
         ELSE
 
86
            FK(I) = 1./(VP(4*I+3)**2)
 
87
         END IF
 
88
   20 CONTINUE
 
89
 
 
90
      VPI    = VP(THREE)
 
91
      VTN(1) = 1.
 
92
 
 
93
      DO 30 K=1,NP
 
94
         VFU = 0.
 
95
         DO 31 L = 1,NC
 
96
            N    = L*3-2
 
97
            N1   = L*4
 
98
            DIX  = IVX(K)-VP(N1+1)
 
99
            DIY  = IVY(K)-VP(N1+2)
 
100
            EPES = DIX*DIX+DIY*DIY
 
101
 
 
102
            IF (BETA.LE.0) THEN
 
103
               EP = EXP(EPES*FK(L))
 
104
               CO = -VP(N1)*EP*2.*FK(L)
 
105
            ELSE
 
106
               EPE = 1.+EPES*FK(L)
 
107
               EP  = EPE**(-BETA)
 
108
               EP2 = EPE**(-BETA-1.)
 
109
               CO  = VP(N1)*BETA*EP2*2*FK(L)
 
110
            END IF
 
111
 
 
112
            VTN(N)   = EP
 
113
            VTN(N+1) = CO*DIX
 
114
            VTN(N+2) = CO*DIY
 
115
            VFU      = VP(N1)*EP+VFU
 
116
   31    CONTINUE
 
117
 
 
118
         VFU = VFU+VPI
 
119
         DO 32 I = 1,NIN
 
120
            VTP    = VTN(I)*WEI(K)
 
121
            VC(I)  = VC(I)+(VZ(K)-VFU)*VTP
 
122
            DO 36 J = 1,I
 
123
               RMT(I,J) = RMT(I,J)+VTN(J)*VTP
 
124
   36       CONTINUE
 
125
   32    CONTINUE
 
126
   30 CONTINUE
 
127
 
 
128
      DO 40 I = 2,NIN
 
129
         DO 41 J = 1,I-1
 
130
            RMT(J,I) = RMT(I,J)
 
131
   41    CONTINUE
 
132
   40 CONTINUE
 
133
 
 
134
      DO 50 I=1,NIN
 
135
         RMT(I,I) = RMT(I,I)*(1+FS**2)
 
136
   50 CONTINUE
 
137
 
 
138
      NCD = NIND
 
139
      CALL LISIB(RMT,VC,NIN,NCD,SIG)
 
140
 
 
141
      IF (NCD.GT.0) THEN
 
142
         DO 60 I = 1,NC
 
143
             DO 61 J = 4,6
 
144
               INK     = J+4*(I-1)
 
145
               INK1    = J+3*(I-1)
 
146
               VP(INK) = VC(INK1)*VPES(J)+VP(INK)
 
147
               IF(ABS(VP(INK)).GT.1000..AND.J.NE.4) NCD=-1
 
148
   61        CONTINUE
 
149
   60    CONTINUE
 
150
 
 
151
         IF (NCD.GT.0) THEN
 
152
            SQM = 0.
 
153
            DO 63 I = 1,NP
 
154
               VG = VP(THREE)
 
155
               DO 64 J = 1,NC
 
156
                  N   = J*4
 
157
                  EEE = ((VP(N+1)-IVX(I))**2 + (VP(N+2)-IVY(I))**2)/
 
158
     2                                         (VP(N+3)**2)
 
159
                  IF(BETA.LE.0.) THEN
 
160
                     VG = VG+VP(N)*EXP(-EEE*4*ALOG(2.))
 
161
                  ELSE
 
162
                     VG = VG+VP(N)*(1+EEE)**(-BETA)
 
163
                  END IF
 
164
   64          CONTINUE
 
165
               SQM = SQM+(VZ(I)-VG)**2*WEI(I)
 
166
   63       CONTINUE
 
167
            SQM = SQM/(NP-NIN)
 
168
         END IF
 
169
      END IF
 
170
 
 
171
      IF(NCD.LT.1) LG=1
 
172
 
 
173
      RETURN
 
174
      END
 
175
 
 
176
      SUBROUTINE ELMRFV(IVX,IVY,VZ,NP,VP,FS,VPES,NC,BETA,SQM,LG,WEI,
 
177
     2                  SIG)
 
178
C+++
 
179
C.PURPOSE: Moffat o Gauss sigma fisso - piano fisso
 
180
C---
 
181
      IMPLICIT     NONE
 
182
      INTEGER      NST
 
183
      INTEGER      NIND
 
184
      PARAMETER    (NST=40)
 
185
      PARAMETER    (NIND=4*NST+3)
 
186
C
 
187
      INTEGER      IVX(1), IVY(1)
 
188
      REAL         VZ(1)
 
189
      INTEGER      NP
 
190
      REAL         VP(1)
 
191
      REAL         FS
 
192
      REAL         VPES(7)
 
193
      INTEGER      NC
 
194
      REAL         BETA
 
195
      REAL         SQM
 
196
      INTEGER      LG, THREE
 
197
      REAL         WEI(1)
 
198
      REAL         SIG
 
199
C
 
200
      REAL         VC, VTP, VPI
 
201
      REAL         CO
 
202
      INTEGER      N1
 
203
      INTEGER      INK, INK1, NCD
 
204
      REAL         EPE, EP2, EEE
 
205
      REAL         RMT
 
206
      REAL         VTN, VFZ, VFF, VFU, VG
 
207
      REAL         DIX, DIY
 
208
      REAL         EP, EPES
 
209
      INTEGER      NIN
 
210
      INTEGER      I, J, K, L, N
 
211
      REAL         FK(NST)
 
212
C
 
213
      COMMON      /SUFR/RMT(NIND,NIND),VTN(NIND),VFZ(NIND),VFF(NIND),
 
214
     2                  VC(NIND)
 
215
C
 
216
C *** start the code
 
217
      THREE = 3
 
218
      NIN=NC*4
 
219
 
 
220
      DO 10 I=1,NIN
 
221
         VFF(I) = 0
 
222
         VFZ(I) = 0
 
223
         VC(I)  = 0
 
224
         DO 11 J = 1,NIN
 
225
            RMT(I,J) = 0
 
226
   11    CONTINUE
 
227
   10 CONTINUE
 
228
 
 
229
      DO 20 I = 1,NC
 
230
         IF(BETA.LE.0.) THEN
 
231
            FK(I) = -4*ALOG(2.)/(VP(4*I+3)**2)
 
232
         ELSE
 
233
            FK(I) = 1./(VP(4*I+3)**2)
 
234
         END IF
 
235
   20 CONTINUE
 
236
 
 
237
      VPI    = VP(THREE)
 
238
      VTN(1) = 1.
 
239
 
 
240
      DO 30 K=1,NP
 
241
         VFU = 0.
 
242
         DO 31 L = 1,NC
 
243
            N    = L*4-3
 
244
            N1   = L*4
 
245
            DIX  = IVX(K)-VP(N1+1)
 
246
            DIY  = IVY(K)-VP(N1+2)
 
247
            EPES = DIX*DIX+DIY*DIY
 
248
 
 
249
            IF (BETA.LE.0) THEN
 
250
               EP = EXP(EPES*FK(L))
 
251
               CO = -VP(N1)*EP*2.*FK(L)
 
252
            ELSE
 
253
               EPE = 1.+EPES*FK(L)
 
254
               EP  = EPE**(-BETA)
 
255
               EP2 = EPE**(-BETA-1.)
 
256
               CO  = VP(N1)*BETA*EP2*2*FK(L)
 
257
            END IF
 
258
 
 
259
            VTN(N)   = EP
 
260
            VTN(N+1) = CO*DIX
 
261
            VTN(N+2) = CO*DIY
 
262
            VTN(N+3) = CO*EPES/VP(N1+3)
 
263
            VFU      = VP(N1)*EP+VFU
 
264
   31    CONTINUE
 
265
 
 
266
         VFU = VFU+VPI
 
267
         DO 32 I = 1,NIN
 
268
            VTP    = VTN(I)*WEI(K)
 
269
            VC(I)  = VC(I)+(VZ(K)-VFU)*VTP
 
270
            DO 36 J = 1,I
 
271
               RMT(I,J) = RMT(I,J)+VTN(J)*VTP
 
272
   36       CONTINUE
 
273
   32    CONTINUE
 
274
   30 CONTINUE
 
275
 
 
276
      DO 40 I = 2,NIN
 
277
         DO 41 J = 1,I-1
 
278
            RMT(J,I) = RMT(I,J)
 
279
   41    CONTINUE
 
280
   40 CONTINUE
 
281
 
 
282
      DO 50 I=1,NIN
 
283
         RMT(I,I) = RMT(I,I)*(1+FS**2)
 
284
   50 CONTINUE
 
285
 
 
286
      NCD = NIND
 
287
      CALL LISIB(RMT,VC,NIN,NCD,SIG)
 
288
 
 
289
      IF (NCD.GT.0) THEN
 
290
         DO 60 I = 1,NC
 
291
             DO 61 J = 4,7
 
292
               INK     = J+4*(I-1)
 
293
               INK1    = INK-3
 
294
               VP(INK) = VC(INK1)*VPES(J)+VP(INK)
 
295
               IF(ABS(VP(INK)).GT.1000..AND.J.NE.4) NCD=-1
 
296
   61        CONTINUE
 
297
   60    CONTINUE
 
298
 
 
299
         IF (NCD.GT.0) THEN
 
300
            SQM = 0.
 
301
            DO 63 I = 1,NP
 
302
               VG = VP(THREE)
 
303
               DO 64 J = 1,NC
 
304
                  N   = J*4
 
305
                  EEE = ((VP(N+1)-IVX(I))**2 + (VP(N+2)-IVY(I))**2)/
 
306
     2                                         (VP(N+3)**2)
 
307
                  IF(BETA.LE.0.) THEN
 
308
                     VG = VG+VP(N)*EXP(-EEE*4*ALOG(2.))
 
309
                  ELSE
 
310
                     VG = VG+VP(N)*(1+EEE)**(-BETA)
 
311
                  END IF
 
312
   64          CONTINUE
 
313
               SQM = SQM+(VZ(I)-VG)**2*WEI(I)
 
314
   63       CONTINUE
 
315
            SQM = SQM/(NP-NIN)
 
316
         END IF
 
317
      END IF
 
318
 
 
319
      IF(NCD.LT.1) LG=1
 
320
 
 
321
      RETURN
 
322
      END
 
323
 
 
324
      SUBROUTINE ELMRR(IVX,IVY,VZ,NP,VP,FS,VPES,NC,BETA,SQM,LG,WEI,
 
325
     2                  SIG)
 
326
C+++
 
327
C.PURPOSE: Moffat o Gauss sigma fisso - piano orizzontale
 
328
C---
 
329
      IMPLICIT     NONE
 
330
      INTEGER      NST
 
331
      INTEGER      NIND
 
332
      PARAMETER    (NST=40)
 
333
      PARAMETER    (NIND=4*NST+3)
 
334
C
 
335
      INTEGER      IVX(1), IVY(1)
 
336
      REAL         VZ(1)
 
337
      INTEGER      NP
 
338
      REAL         VP(1)
 
339
      REAL         FS
 
340
      REAL         VPES(7)
 
341
      INTEGER      NC
 
342
      REAL         BETA
 
343
      REAL         SQM
 
344
      INTEGER      LG, THREE
 
345
      REAL         WEI(1)
 
346
      REAL         SIG
 
347
C
 
348
      REAL         VC, VTP, VPI
 
349
      REAL         CO
 
350
      INTEGER      N1
 
351
      INTEGER      INK, INK1, NCD
 
352
      REAL         EPE, EP2, EEE
 
353
      REAL         RMT
 
354
      REAL         VTN, VFZ, VFF, VFU, VG
 
355
      REAL         DIX, DIY
 
356
      REAL         EP, EPES
 
357
      INTEGER      NIN
 
358
      INTEGER      I, J, K, L, N
 
359
      REAL         FK(NST)
 
360
C
 
361
      COMMON      /SUFR/RMT(NIND,NIND),VTN(NIND),VFZ(NIND),VFF(NIND),
 
362
     2                  VC(NIND)
 
363
C
 
364
C *** start the code
 
365
      THREE = 3
 
366
      NIN=NC*3+1
 
367
 
 
368
      DO 10 I=1,NIN
 
369
         VFF(I) = 0
 
370
         VFZ(I) = 0
 
371
         VC(I)  = 0
 
372
         DO 11 J = 1,NIN
 
373
            RMT(I,J) = 0
 
374
   11    CONTINUE
 
375
   10 CONTINUE
 
376
 
 
377
      DO 20 I = 1,NC
 
378
         IF(BETA.LE.0.) THEN
 
379
            FK(I) = -4*ALOG(2.)/(VP(4*I+3)**2)
 
380
         ELSE
 
381
            FK(I) = 1./(VP(4*I+3)**2)
 
382
         END IF
 
383
   20 CONTINUE
 
384
 
 
385
      VPI    = VP(THREE)
 
386
      VTN(1) = 1.
 
387
 
 
388
      DO 30 K=1,NP
 
389
         VFU = 0.
 
390
         DO 31 L = 1,NC
 
391
            N    = L*3-1
 
392
            N1   = L*4
 
393
            DIX  = IVX(K)-VP(N1+1)
 
394
            DIY  = IVY(K)-VP(N1+2)
 
395
            EPES = DIX*DIX+DIY*DIY
 
396
 
 
397
            IF (BETA.LE.0) THEN
 
398
               EP = EXP(EPES*FK(L))
 
399
               CO = -VP(N1)*EP*2.*FK(L)
 
400
            ELSE
 
401
               EPE = 1.+EPES*FK(L)
 
402
               EP  = EPE**(-BETA)
 
403
               EP2 = EPE**(-BETA-1.)
 
404
               CO  = VP(N1)*BETA*EP2*2*FK(L)
 
405
            END IF
 
406
 
 
407
            VTN(N)   = EP
 
408
            VTN(N+1) = CO*DIX
 
409
            VTN(N+2) = CO*DIY
 
410
            VFU      = VP(N1)*EP+VFU
 
411
   31    CONTINUE
 
412
 
 
413
         VFU = VFU+VPI
 
414
         DO 32 I = 1,NIN
 
415
            VTP    = VTN(I)*WEI(K)
 
416
            VC(I)  = VC(I)+(VZ(K)-VFU)*VTP
 
417
            DO 36 J = 1,I
 
418
               RMT(I,J) = RMT(I,J)+VTN(J)*VTP
 
419
   36       CONTINUE
 
420
   32    CONTINUE
 
421
   30 CONTINUE
 
422
 
 
423
      DO 40 I = 2,NIN
 
424
         DO 41 J = 1,I-1
 
425
            RMT(J,I) = RMT(I,J)
 
426
   41    CONTINUE
 
427
   40 CONTINUE
 
428
 
 
429
      DO 50 I=1,NIN
 
430
         RMT(I,I) = RMT(I,I)*(1+FS**2)
 
431
   50 CONTINUE
 
432
 
 
433
      NCD = NIND
 
434
      CALL LISIB(RMT,VC,NIN,NCD,SIG)
 
435
 
 
436
      IF (NCD.GT.0) THEN
 
437
         VP(THREE) = VC(1)*VPES(3)+VP(THREE)
 
438
         DO 60 I = 1,NC
 
439
             DO 61 J = 4,6
 
440
               INK     = J+4*(I-1)
 
441
               INK1    = J+3*(I-1)-2
 
442
               VP(INK) = VC(INK1)*VPES(J)+VP(INK)
 
443
               IF(ABS(VP(INK)).GT.1000..AND.J.NE.4) NCD=-1
 
444
   61        CONTINUE
 
445
   60    CONTINUE
 
446
 
 
447
         IF (NCD.GT.0) THEN
 
448
            SQM = 0.
 
449
            DO 63 I = 1,NP
 
450
               VG = VP(THREE)
 
451
               DO 64 J = 1,NC
 
452
                  N   = J*4
 
453
                  EEE = ((VP(N+1)-IVX(I))**2 + (VP(N+2)-IVY(I))**2)/
 
454
     2                                         (VP(N+3)**2)
 
455
                  IF(BETA.LE.0.) THEN
 
456
                     VG = VG+VP(N)*EXP(-EEE*4*ALOG(2.))
 
457
                  ELSE
 
458
                     VG = VG+VP(N)*(1+EEE)**(-BETA)
 
459
                  END IF
 
460
   64          CONTINUE
 
461
               SQM = SQM+(VZ(I)-VG)**2*WEI(I)
 
462
   63       CONTINUE
 
463
            SQM = SQM/(NP-NIN)
 
464
         END IF
 
465
      END IF
 
466
 
 
467
      IF(NCD.LT.1) LG=1
 
468
 
 
469
      RETURN
 
470
      END
 
471
 
 
472
      SUBROUTINE ELMRRV(IVX,IVY,VZ,NP,VP,FS,VPES,NC,BETA,SQM,LG,WEI,
 
473
     2                  SIG)
 
474
C
 
475
C     MOFFAT O GAUSS SIGMA VAR. - PIANO ORIZZONTALE
 
476
C
 
477
C-------------
 
478
      IMPLICIT     NONE
 
479
      INTEGER      NST
 
480
      INTEGER      NIND
 
481
      PARAMETER (NST=40)
 
482
      PARAMETER (NIND=NST*4+3)
 
483
C
 
484
      INTEGER      IVX(1), IVY(1)
 
485
      REAL         VZ(1)
 
486
      INTEGER      NP
 
487
      REAL         VP(1)
 
488
      REAL         FS
 
489
      REAL         VPES(7)
 
490
      INTEGER      NC
 
491
      REAL         BETA
 
492
      REAL         SQM
 
493
      INTEGER      LG, THREE
 
494
      REAL         WEI(1)
 
495
      REAL         SIG
 
496
C
 
497
      REAL         VC, VTP, VPI
 
498
      REAL         CO
 
499
      INTEGER      N1
 
500
      INTEGER      INK, INK1, NCD
 
501
      REAL         EPE, EP2, EEE
 
502
      REAL         RMT
 
503
      REAL         VTN, VFZ, VFF, VFU, VG
 
504
      REAL         DIX, DIY
 
505
      REAL         EP, EPES
 
506
      INTEGER      NIN
 
507
      INTEGER      I, J, K, L, N
 
508
      REAL         FK(NST)
 
509
C
 
510
      COMMON /SUFR/RMT(NIND,NIND),VTN(NIND),VFZ(NIND),VFF(NIND),
 
511
     2             VC(NIND)
 
512
 
 
513
      THREE = 3
 
514
      NIN = NC*4+1
 
515
      DO 10 I=1,NIN
 
516
         VFF(I)=0
 
517
         VFZ(I)=0
 
518
         VC(I)=0
 
519
         DO 11 J=1,NIN
 
520
            RMT(I,J)=0
 
521
   11    CONTINUE
 
522
   10 CONTINUE
 
523
C
 
524
      DO 20 I=1,NC
 
525
         IF (BETA.LE.0.) THEN
 
526
            FK(I) = -4*ALOG(2.)/(VP(4*I+3)**2)
 
527
         ELSE
 
528
            FK(I) = 1./(VP(4*I+3)**2)
 
529
         END IF
 
530
   20 CONTINUE
 
531
C
 
532
      VPI    = VP(THREE)
 
533
      VTN(1) = 1.
 
534
      DO 30 K=1,NP
 
535
         VFU = 0.
 
536
         DO 31 L=1,NC
 
537
            N    = L*4-2
 
538
            N1   = L*4
 
539
            DIX  = IVX(K)-VP(N1+1)
 
540
            DIY  = IVY(K)-VP(N1+2)
 
541
            EPES = DIX*DIX+DIY*DIY
 
542
            IF (BETA.LE.0) THEN
 
543
               EP  = EXP(EPES*FK(L))
 
544
               CO  = -VP(N1)*EP*2.*FK(L)
 
545
            ELSE
 
546
               EPE = 1.+EPES*FK(L)
 
547
               EP  = EPE**(-BETA)
 
548
               EP2 = EPE**(-BETA-1.)
 
549
               CO  = VP(N1)*BETA*EP2*2*FK(L)
 
550
            END IF
 
551
C
 
552
            VTN(N)   = EP
 
553
            VTN(N+1) = CO*DIX
 
554
            VTN(N+2) = CO*DIY 
 
555
            VTN(N+3) = CO*EPES/VP(N1+3)
 
556
            VFU      = VP(N1)*EP+VFU
 
557
   31    CONTINUE
 
558
C
 
559
         VFU = VFU+VPI
 
560
         DO 32 I=1,NIN
 
561
            VTP   = VTN(I)*WEI(K)
 
562
            VC(I) = VC(I)+(VZ(K)-VFU)*VTP
 
563
            DO 132 J=1,I
 
564
               RMT(I,J) = RMT(I,J)+VTN(J)*VTP
 
565
  132       CONTINUE
 
566
   32    CONTINUE
 
567
   30 CONTINUE
 
568
C 
 
569
      DO 40 I=2,NIN
 
570
         DO 41 J=1,I-1
 
571
            RMT(J,I) = RMT(I,J)
 
572
   41    CONTINUE
 
573
   40 CONTINUE
 
574
C
 
575
      DO 50 I=1,NIN
 
576
         RMT(I,I) = RMT(I,I)*(1+FS**2)
 
577
   50 CONTINUE
 
578
C
 
579
      NCD=NIND
 
580
      CALL LISIB(RMT,VC,NIN,NCD,SIG)
 
581
      IF (NCD.GT.0) THEN
 
582
         VP(THREE) = VC(1)*VPES(3)+VP(THREE)
 
583
         DO 61 I=1,NC
 
584
            DO 161 J=4,7
 
585
               INK     = J+4*(I-1)
 
586
               INK1    = INK-2
 
587
               VP(INK) = VC(INK1)*VPES(J)+VP(INK)
 
588
               IF (ABS(VP(INK)).GT.1000..AND.J.NE.4) NCD=-1
 
589
  161       CONTINUE
 
590
   61    CONTINUE
 
591
C
 
592
         IF (NCD.GT.0) THEN 
 
593
            SQM = 0.
 
594
            DO 62 I=1,NP
 
595
               VG=VP(THREE)
 
596
               DO 162 J=1,NC
 
597
                  N   = J*4
 
598
                  EEE = ((VP(N+1)-IVX(I))**2+(VP(N+2)-IVY(I))**2)/
 
599
     2                   (VP(N+3)**2)
 
600
                  IF (BETA.LE.0.) THEN
 
601
                     VG = VG+VP(N)*EXP(-EEE*4*ALOG(2.))
 
602
                  ELSE
 
603
                     VG = VG+VP(N)*(1+EEE)**(-BETA)
 
604
                  END IF
 
605
  162          CONTINUE
 
606
               SQM = SQM+(VZ(I)-VG)**2*WEI(I)
 
607
   62       CONTINUE
 
608
            SQM = SQM/(NP-NIN)
 
609
         END IF
 
610
      END IF
 
611
C
 
612
      IF (NCD.LT.1) LG=1
 
613
      RETURN
 
614
      END
 
615
 
 
616
      SUBROUTINE ELMRPF(IVX,IVY,VZ,NP,VP,FS,VPES,NC,BETA,SQM,LG,WEI,SIG)
 
617
C+++
 
618
C.PURPOSE: Moffat o Gauss sigma fisso - piano orizzontale - posiz. fissa
 
619
C---  
 
620
      IMPLICIT     NONE
 
621
      INTEGER      NST
 
622
      INTEGER      NIND
 
623
      PARAMETER (NST=40)
 
624
      PARAMETER (NIND=4*NST+3)
 
625
C
 
626
      INTEGER      IVX(1), IVY(1)
 
627
      REAL         VZ(1)
 
628
      INTEGER      NP
 
629
      REAL         VP(1)
 
630
      REAL         FS
 
631
      REAL         VPES(7)
 
632
      INTEGER      NC
 
633
      REAL         BETA
 
634
      REAL         SQM
 
635
      INTEGER      LG, THREE
 
636
      REAL         WEI(1)
 
637
      REAL         SIG
 
638
C
 
639
      REAL         VC, VTP
 
640
      INTEGER      N1
 
641
      INTEGER      NCD
 
642
      REAL         EPE
 
643
      REAL         RMT
 
644
      REAL         VTN, VFZ, VFF, VFU, VG
 
645
      REAL         DIX, DIY
 
646
      REAL         EP, EPES,EEE
 
647
      INTEGER      NIN
 
648
      INTEGER      I, J, K, L, N
 
649
      REAL         FK(NST)
 
650
C
 
651
      COMMON    /SUFR/RMT(NIND,NIND),VTN(NIND),VFZ(NIND),
 
652
     2                VFF(NIND),VC(NIND)
 
653
C
 
654
      THREE = 3
 
655
      NIN = NC+1
 
656
      DO 10 I = 1,NIN
 
657
         VC(I)=0
 
658
         DO 20  J = 1,NIN
 
659
            RMT(I,J)=0
 
660
   20    CONTINUE
 
661
   10 CONTINUE
 
662
C
 
663
      DO 30 I = 1,NC
 
664
         IF (BETA.LE.0.) THEN
 
665
            FK(I) = -4*ALOG(2.)/(VP(4*I+3)**2)
 
666
         ELSE
 
667
            FK(I) = 1./(VP(4*I+3)**2)
 
668
         END IF
 
669
   30 CONTINUE
 
670
C
 
671
      VTN(1) = 1.
 
672
      DO 40  K = 1,NP
 
673
         VFU = 0.
 
674
         DO 41 L=1,NC
 
675
            N1   = L*4
 
676
            DIX  = IVX(K)-VP(N1+1)
 
677
            DIY  = IVY(K)-VP(N1+2)
 
678
            EPES = DIX*DIX+DIY*DIY
 
679
            IF (BETA.LE.0) THEN
 
680
               EP  = EXP(EPES*FK(L))
 
681
            ELSE
 
682
               EPE = 1.+EPES*FK(L)
 
683
               EP  = EPE**(-BETA)
 
684
            END IF
 
685
            VTN(L+1) = EP
 
686
   41    CONTINUE
 
687
C
 
688
         DO 42 I = 1,NIN
 
689
            VTP   = VTN(I)*WEI(K)
 
690
            VC(I) = VC(I)+VZ(K)*VTP
 
691
            DO 43 J=1,I
 
692
               RMT(I,J) = RMT(I,J)+VTN(J)*VTP
 
693
   43       CONTINUE
 
694
   42    CONTINUE
 
695
   40 CONTINUE
 
696
C
 
697
      DO 50 I = 2,NIN
 
698
         DO 51 J=1,I-1
 
699
            RMT(J,I)=RMT(I,J)
 
700
   51    CONTINUE
 
701
   50 CONTINUE
 
702
C
 
703
      NCD = NIND
 
704
      CALL LISIB(RMT,VC,NIN,NCD,SIG)
 
705
      IF (NCD.GT.0) THEN
 
706
         VP(THREE) = VC(1)
 
707
         DO 60 I=1,NC
 
708
            VP(I*4) = VC(I+1)
 
709
   60    CONTINUE
 
710
         SQM=0.
 
711
         DO 70 I=1,NP
 
712
            VG = VP(THREE)
 
713
            DO 71 J=1,NC
 
714
               N   = J*4
 
715
               EEE = ((VP(N+1)-IVX(I))**2+(VP(N+2)-IVY(I))**2)/
 
716
     2                (VP(N+3)**2)
 
717
               IF (BETA.LE.0.) THEN
 
718
                  VG = VG+VP(N)*EXP(-EEE*4*ALOG(2.))
 
719
               ELSE
 
720
                  VG = VG+VP(N)*(1+EEE)**(-BETA)
 
721
               END IF
 
722
   71       CONTINUE
 
723
            SQM=SQM+(VZ(I)-VG)**2*WEI(I)
 
724
   70    CONTINUE
 
725
         SQM=SQM/(NP-NIN)
 
726
      END IF
 
727
C
 
728
      IF (NCD.LT.1) LG=1
 
729
      RETURN
 
730
      END
 
731
 
 
732
      SUBROUTINE ELMRPV(IVX,IVY,VZ,NP,VP,FS,VPES,NC,BETA,SQM,LG,WEI,
 
733
     2                  SIG)
 
734
C+++
 
735
C.PURPOSE: Moffat o Gauss sigma fisso - piano orizzontale - posiz. fissa
 
736
C---
 
737
      IMPLICIT     NONE
 
738
      INTEGER      NST
 
739
      INTEGER      NIND
 
740
      PARAMETER    (NST=40)
 
741
      PARAMETER    (NIND=4*NST+3)
 
742
C
 
743
      INTEGER      IVX(1), IVY(1)
 
744
      REAL         VZ(1)
 
745
      INTEGER      NP
 
746
      REAL         VP(1)
 
747
      REAL         FS
 
748
      REAL         VPES(7)
 
749
      INTEGER      NC
 
750
      REAL         BETA
 
751
      REAL         SQM
 
752
      INTEGER      LG, THREE
 
753
      REAL         WEI(1)
 
754
      REAL         SIG
 
755
C
 
756
      REAL         VC, VTP, VPI
 
757
      REAL         CO
 
758
      INTEGER      N1
 
759
      INTEGER      INK, INK1, NCD
 
760
      REAL         EPE, EP2, EEE
 
761
      REAL         RMT
 
762
      REAL         VTN, VFZ, VFF, VFU, VG
 
763
      REAL         DIX, DIY
 
764
      REAL         EP, EPES
 
765
      INTEGER      NIN
 
766
      INTEGER      I, J, K, L, N
 
767
      REAL         FK(NST)
 
768
C
 
769
      COMMON      /SUFR/RMT(NIND,NIND),VTN(NIND),VFZ(NIND),VFF(NIND),
 
770
     2                  VC(NIND)
 
771
C
 
772
C *** start the code
 
773
      THREE = 3 
 
774
      NIN=NC*2+1
 
775
 
 
776
      DO 10 I=1,NIN
 
777
         VFF(I) = 0
 
778
         VFZ(I) = 0
 
779
         VC(I)  = 0
 
780
         DO 11 J = 1,NIN
 
781
            RMT(I,J) = 0
 
782
   11    CONTINUE
 
783
   10 CONTINUE
 
784
 
 
785
      DO 20 I = 1,NC
 
786
         IF(BETA.LE.0.) THEN
 
787
            FK(I) = -4*ALOG(2.)/(VP(4*I+3)**2)
 
788
         ELSE
 
789
            FK(I) = 1./(VP(4*I+3)**2)
 
790
         END IF
 
791
   20 CONTINUE
 
792
 
 
793
      VPI    = VP(THREE)
 
794
      VTN(1) = 1.
 
795
 
 
796
      DO 30 K=1,NP
 
797
         VFU = 0.
 
798
         DO 31 L = 1,NC
 
799
            N    = L*2
 
800
            N1   = L*4
 
801
            DIX  = IVX(K)-VP(N1+1)
 
802
            DIY  = IVY(K)-VP(N1+2)
 
803
            EPES = DIX*DIX+DIY*DIY
 
804
 
 
805
            IF (BETA.LE.0) THEN
 
806
               EP = EXP(EPES*FK(L))
 
807
               CO = -VP(N1)*EP*2.*FK(L)
 
808
            ELSE
 
809
               EPE = 1.+EPES*FK(L)
 
810
               EP  = EPE**(-BETA)
 
811
               EP2 = EPE**(-BETA-1.)
 
812
               CO  = VP(N1)*BETA*EP2*2*FK(L)
 
813
            END IF
 
814
 
 
815
            VTN(N)   = EP
 
816
            VTN(N+1) = CO*EPES/VP(N1+3)
 
817
            VFU      = VP(N1)*EP+VFU
 
818
   31    CONTINUE
 
819
 
 
820
         VFU = VFU+VPI
 
821
         DO 32 I = 1,NIN
 
822
            VTP    = VTN(I)*WEI(K)
 
823
            VC(I)  = VC(I)+(VZ(K)-VFU)*VTP
 
824
            DO 36 J = 1,I
 
825
               RMT(I,J) = RMT(I,J)+VTN(J)*VTP
 
826
   36       CONTINUE
 
827
   32    CONTINUE
 
828
   30 CONTINUE
 
829
 
 
830
      DO 40 I = 2,NIN
 
831
         DO 41 J = 1,I-1
 
832
            RMT(J,I) = RMT(I,J)
 
833
   41    CONTINUE
 
834
   40 CONTINUE
 
835
 
 
836
      DO 50 I=1,NIN
 
837
         RMT(I,I) = RMT(I,I)*(1+FS**2)
 
838
   50 CONTINUE
 
839
 
 
840
      NCD = NIND
 
841
      CALL LISIB(RMT,VC,NIN,NCD,SIG)
 
842
 
 
843
      IF (NCD.GT.0) THEN
 
844
         VP(THREE) = VC(1)*VPES(3)+VP(THREE)
 
845
         DO 60 I = 1,NC
 
846
             DO 61 J = 4,7,3
 
847
               INK     = J+4*(I-1)
 
848
               INK1    = INK/2
 
849
               VP(INK) = VC(INK1)*VPES(J)+VP(INK)
 
850
               IF(ABS(VP(INK)).GT.1000..AND.J.NE.4) NCD=-1
 
851
   61        CONTINUE
 
852
   60    CONTINUE
 
853
 
 
854
         IF (NCD.GT.0) THEN
 
855
            SQM = 0.
 
856
            DO 63 I = 1,NP
 
857
               VG = VP(THREE)
 
858
               DO 64 J = 1,NC
 
859
                  N   = J*4
 
860
                  EEE = ((VP(N+1)-IVX(I))**2 + (VP(N+2)-IVY(I))**2)/
 
861
     2                                         (VP(N+3)**2)
 
862
                  IF(BETA.LE.0.) THEN
 
863
                     VG = VG+VP(N)*EXP(-EEE*4*ALOG(2.))
 
864
                  ELSE
 
865
                     VG = VG+VP(N)*(1+EEE)**(-BETA)
 
866
                  END IF
 
867
   64          CONTINUE
 
868
               SQM = SQM+(VZ(I)-VG)**2*WEI(I)
 
869
   63       CONTINUE
 
870
            SQM = SQM/(NP-NIN)
 
871
         END IF
 
872
      END IF
 
873
 
 
874
      IF(NCD.LT.1) LG=1
 
875
 
 
876
      RETURN
 
877
      END
 
878
 
 
879
      SUBROUTINE ELMRV(IVX,IVY,VZ,NP,VP,FS,VPES,NC,BETA,SQM,LG,WEI,SIG)
 
880
C+++
 
881
C.PURPOSE: Moffat o Gauss sigma var. - piano inclinato
 
882
C---
 
883
      IMPLICIT     NONE
 
884
      INTEGER      NST
 
885
      INTEGER      NIND
 
886
      PARAMETER (NST=40)
 
887
      PARAMETER (NIND=4*NST+3)
 
888
C
 
889
      INTEGER      IVX(1), IVY(1)
 
890
      REAL         VZ(1)
 
891
      INTEGER      NP
 
892
      REAL         VP(1)
 
893
      REAL         FS
 
894
      REAL         VPES(7)
 
895
      INTEGER      NC
 
896
      REAL         BETA
 
897
      REAL         SQM
 
898
      INTEGER      LG, THREE, TWO
 
899
      REAL         WEI(1)
 
900
      REAL         SIG
 
901
C
 
902
      REAL         VC, VTP, VPI
 
903
      REAL         CO
 
904
      INTEGER      N1
 
905
      INTEGER      INK, NCD
 
906
      REAL         EPE, EP2, EEE
 
907
      REAL         RMT
 
908
      REAL         VTN, VFZ, VFF, VFU, VG
 
909
      REAL         DIX, DIY
 
910
      REAL         EP, EPES
 
911
      INTEGER      NIN
 
912
      INTEGER      I, J, K, L, N
 
913
      REAL         FK(NST)
 
914
C
 
915
      COMMON    /SUFR/RMT(NIND,NIND),VTN(NIND),VFZ(NIND),
 
916
     2                 VFF(NIND),VC(NIND)
 
917
C
 
918
      TWO   = 2
 
919
      THREE = 3
 
920
      NIN = NC*4+3
 
921
      DO 10 I = 1,NIN
 
922
         VFF(I) = 0
 
923
         VFZ(I) = 0
 
924
         VC(I)  = 0
 
925
         DO 11 J = 1,NIN
 
926
            RMT(I,J) = 0
 
927
   11    CONTINUE
 
928
   10 CONTINUE
 
929
C
 
930
      DO 20 I = 1,NC
 
931
         IF (BETA.LE.0.) THEN
 
932
            FK(I) = -4*ALOG(2.)/(VP(4*I+3)**2)
 
933
         ELSE
 
934
            FK(I) = 1./(VP(4*I+3)**2)
 
935
         END IF
 
936
   20 CONTINUE
 
937
C
 
938
      VTN(THREE) = 1.
 
939
      DO 30 K = 1,NP
 
940
         VTN(1) = IVX(K)
 
941
         VTN(2) = IVY(K)
 
942
         VPI    = VP(1)*IVX(K)+VP(TWO)*IVY(K)+VP(THREE)
 
943
         VFU    = 0.
 
944
         DO 31 L = 1,NC
 
945
            N1   = L*4
 
946
            DIX  = IVX(K)-VP(N1+1)
 
947
            DIY  = IVY(K)-VP(N1+2)
 
948
            EPES = DIX*DIX+DIY*DIY
 
949
            IF (BETA.LE.0) THEN
 
950
               EP = EXP(EPES*FK(L))
 
951
               CO = -VP(N1)*EP*2.*FK(L)
 
952
            ELSE
 
953
               EPE = 1.+EPES*FK(L)
 
954
               EP  = EPE**(-BETA)
 
955
               EP2 = EPE**(-BETA-1.)
 
956
               CO  = VP(N1)*BETA*EP2*2*FK(L)
 
957
            END IF
 
958
            VTN(N1)   = EP
 
959
            VTN(N1+1) = CO*DIX
 
960
            VTN(N1+2) = CO*DIY
 
961
            VTN(N1+3) = CO*EPES/VP(N1+3)
 
962
            VFU       = VP(N1)*EP+VFU
 
963
   31    CONTINUE
 
964
C
 
965
         VFU = VFU+VPI
 
966
         DO 32 I = 1,NIN
 
967
            VTP   = VTN(I)*WEI(K)
 
968
            VC(I) = VC(I)+(VZ(K)-VFU)*VTP
 
969
            DO 33 J = 1,I
 
970
               RMT(I,J) = RMT(I,J)+VTN(J)*VTP
 
971
   33       CONTINUE
 
972
   32     CONTINUE
 
973
   30 CONTINUE
 
974
C
 
975
      DO 40 I = 2,NIN
 
976
         DO 41 J = 1,I-1
 
977
            RMT(J,I) = RMT(I,J)
 
978
   41    CONTINUE
 
979
   40 CONTINUE
 
980
 
 
981
      DO 50 I = 1,NIN
 
982
         RMT(I,I) = RMT(I,I)*(1+FS**2)
 
983
   50 CONTINUE
 
984
C
 
985
      NCD = NIND
 
986
      CALL LISIB(RMT,VC,NIN,NCD,SIG)
 
987
      IF (NCD.GT.0) THEN
 
988
         DO 60 I = 1,3
 
989
            VP(I) = VC(I)*VPES(I)+VP(I)
 
990
   60    CONTINUE
 
991
C
 
992
         DO 70 I = 1,NC
 
993
            DO 71 J = 4,7
 
994
              INK     = J+4*(I-1)
 
995
               VP(INK) = VC(INK)*VPES(J)+VP(INK)
 
996
               IF(ABS(VP(INK)).GT.1000..AND.J.NE.4) NCD=-1
 
997
   71       CONTINUE
 
998
   70    CONTINUE
 
999
C
 
1000
         IF (NCD.GT.0) THEN
 
1001
            SQM = 0.
 
1002
            DO 80 I = 1,NP
 
1003
               VG = VP(1)*IVX(I)+VP(TWO)*IVY(I)+VP(THREE)
 
1004
               DO 81 J = 1,NC
 
1005
                  N   = J*4
 
1006
                  EEE = ((VP(N+1)-IVX(I))**2+(VP(N+2)-IVY(I))**2)/
 
1007
     2                   (VP(N+3)**2)
 
1008
                  IF (BETA.LE.0.) THEN
 
1009
                     VG = VG+VP(N)*EXP(-EEE*4*ALOG(2.))
 
1010
                  ELSE
 
1011
                     VG = VG+VP(N)*(1+EEE)**(-BETA)
 
1012
                  END IF
 
1013
   81          CONTINUE
 
1014
               SQM = SQM+(VZ(I)-VG)**2*WEI(I)
 
1015
   80       CONTINUE
 
1016
            SQM = SQM/(NP-NIN)
 
1017
         END IF
 
1018
      END IF
 
1019
C
 
1020
      IF (NCD.LT.1) LG=1
 
1021
      RETURN
 
1022
      END
 
1023
 
 
1024
      SUBROUTINE ELMRH(IVX,IVY,VZ,NP,VP,FS,VPES,NC,BETA,SQM,LG,WEI,SIG)
 
1025
C+++
 
1026
C.PURPOSE: Moffat o Gauss sigma fisso - piano fisso - posiz. fissa
 
1027
C---
 
1028
      IMPLICIT     NONE
 
1029
      INTEGER      NST
 
1030
      INTEGER      NIND
 
1031
      PARAMETER  (NST=40)
 
1032
      PARAMETER  (NIND=4*NST+3)
 
1033
C
 
1034
      INTEGER      IVX(1), IVY(1)
 
1035
      REAL         VZ(1)
 
1036
      INTEGER      NP
 
1037
      REAL         VP(1)
 
1038
      REAL         FS
 
1039
      REAL         VPES(7)
 
1040
      INTEGER      NC
 
1041
      REAL         BETA
 
1042
      REAL         SQM
 
1043
      INTEGER      LG, THREE
 
1044
      REAL         WEI(1)
 
1045
      REAL         SIG
 
1046
C
 
1047
      REAL         VC, VTP
 
1048
      INTEGER      N1
 
1049
      INTEGER      NCD
 
1050
      REAL         EPE, EEE
 
1051
      REAL         RMT
 
1052
      REAL         VTN, VFZ, VFF, VFU, VG
 
1053
      REAL         DIX, DIY
 
1054
      REAL         EP, EPES
 
1055
      INTEGER      NIN
 
1056
      INTEGER      I, J, K, L, N
 
1057
      REAL         FK(NST)
 
1058
C
 
1059
      COMMON    /SUFR/RMT(NIND,NIND),VTN(NIND),VFZ(NIND),
 
1060
     2                VFF(NIND),VC(NIND)
 
1061
C
 
1062
      THREE = 3
 
1063
      NIN = NC
 
1064
      DO 10 I = 1,NIN
 
1065
         VC(I) = 0
 
1066
         DO 11 J = 1,NIN
 
1067
            RMT(I,J) = 0
 
1068
   11    CONTINUE
 
1069
   10 CONTINUE
 
1070
C
 
1071
      DO 20 I = 1,NC
 
1072
         IF (BETA.LE.0.) THEN
 
1073
            FK(I) = -4*ALOG(2.)/(VP(4*I+3)**2)
 
1074
         ELSE
 
1075
            FK(I) = 1./(VP(4*I+3)**2)
 
1076
         END IF
 
1077
   20 CONTINUE
 
1078
C
 
1079
      DO 30 K = 1,NP
 
1080
         VFU = 0.
 
1081
         DO 31 L = 1,NC
 
1082
            N1   = L*4
 
1083
            DIX  = IVX(K)-VP(N1+1)
 
1084
            DIY  = IVY(K)-VP(N1+2)
 
1085
            EPES = DIX*DIX+DIY*DIY
 
1086
            IF (BETA.LE.0) THEN
 
1087
               EP = EXP(EPES*FK(L))
 
1088
            ELSE
 
1089
               EPE = 1.+EPES*FK(L)
 
1090
               EP  = EPE**(-BETA)
 
1091
            END IF
 
1092
            VTN(L) = EP
 
1093
   31    CONTINUE
 
1094
C
 
1095
         DO 32 I = 1,NIN
 
1096
            VTP   = VTN(I)*WEI(K)
 
1097
            VC(I) = VC(I)+VZ(K)*VTP
 
1098
            DO 33 J = 1,I
 
1099
               RMT(I,J) = RMT(I,J)+VTN(J)*VTP
 
1100
   33       CONTINUE
 
1101
   32    CONTINUE
 
1102
   30 CONTINUE
 
1103
C
 
1104
      DO 40 I = 2,NIN
 
1105
         DO 41 J = 1,I-1
 
1106
            RMT(J,I) = RMT(I,J)
 
1107
   41    CONTINUE
 
1108
   40 CONTINUE
 
1109
C
 
1110
      NCD = NIND
 
1111
      CALL LISIB(RMT,VC,NIN,NCD,SIG)
 
1112
      IF (NCD.GT.0) THEN
 
1113
         DO 50 I = 1,NC
 
1114
            VP(I*4) = VC(I)
 
1115
   50    CONTINUE
 
1116
C
 
1117
         SQM = 0.
 
1118
         DO 60 I = 1,NP
 
1119
            VG = VP(THREE)   
 
1120
            DO 61 J = 1,NC
 
1121
               N   = J*4  
 
1122
               EEE = ((VP(N+1)-IVX(I))**2+(VP(N+2)-IVY(I))**2)/
 
1123
     2                (VP(N+3)**2)
 
1124
               IF (BETA.LE.0.) THEN
 
1125
                  VG = VG+VP(N)*EXP(-EEE*4*ALOG(2.))
 
1126
               ELSE
 
1127
                  VG = VG+VP(N)*(1+EEE)**(-BETA)
 
1128
               END IF
 
1129
   61       CONTINUE
 
1130
            SQM = SQM+(VZ(I)-VG)**2*WEI(I)
 
1131
   60    CONTINUE
 
1132
         SQM = SQM/(NP-NIN)
 
1133
      END IF
 
1134
C
 
1135
      IF (NCD.LT.1) LG=1
 
1136
      RETURN
 
1137
      END
 
1138
 
 
1139
      SUBROUTINE ELMRX(IVX,IVY,VZ,NP,VP,FS,VPES,NC,BETA,SQM,LG,WEI,
 
1140
     2                  SIG)
 
1141
C+++
 
1142
C.PURPOSE: Moffat o Gauss sigma fisso - piano orizzontale
 
1143
C---
 
1144
      IMPLICIT     NONE
 
1145
      INTEGER      NST
 
1146
      INTEGER      NIND
 
1147
      PARAMETER    (NST=40)
 
1148
      PARAMETER    (NIND=4*NST+3)
 
1149
C
 
1150
      INTEGER      IVX(1), IVY(1)
 
1151
      REAL         VZ(1)
 
1152
      INTEGER      NP
 
1153
      REAL         VP(1)
 
1154
      REAL         FS
 
1155
      REAL         VPES(7)
 
1156
      INTEGER      NC
 
1157
      REAL         BETA
 
1158
      REAL         SQM
 
1159
      INTEGER      LG, THREE, TWO
 
1160
      REAL         WEI(1)
 
1161
      REAL         SIG
 
1162
C
 
1163
      REAL         VC, VTP, VPI
 
1164
      REAL         CO
 
1165
      INTEGER      N1
 
1166
      INTEGER      INK, INK1, NCD
 
1167
      REAL         EPE, EP2, EEE
 
1168
      REAL         RMT
 
1169
      REAL         VTN, VFZ, VFF, VFU, VG
 
1170
      REAL         DIX, DIY
 
1171
      REAL         EP, EPES
 
1172
      INTEGER      NIN
 
1173
      INTEGER      I, J, K, L, N
 
1174
      REAL         FK(NST)
 
1175
C
 
1176
      COMMON      /SUFR/RMT(NIND,NIND),VTN(NIND),VFZ(NIND),VFF(NIND),
 
1177
     2                  VC(NIND)
 
1178
C
 
1179
C *** start the code
 
1180
      TWO   = 2
 
1181
      THREE = 3
 
1182
      NIN=NC*3+3
 
1183
 
 
1184
      DO 10 I=1,NIN
 
1185
         VFF(I) = 0
 
1186
         VFZ(I) = 0
 
1187
         VC(I)  = 0
 
1188
         DO 11 J = 1,NIN
 
1189
            RMT(I,J) = 0
 
1190
   11    CONTINUE
 
1191
   10 CONTINUE
 
1192
 
 
1193
      DO 20 I = 1,NC
 
1194
         IF(BETA.LE.0.) THEN
 
1195
            FK(I) = -4*ALOG(2.)/(VP(4*I+3)**2)
 
1196
         ELSE
 
1197
            FK(I) = 1./(VP(4*I+3)**2)
 
1198
         END IF
 
1199
   20 CONTINUE
 
1200
 
 
1201
      VTN(THREE) = 1.
 
1202
 
 
1203
      DO 30 K=1,NP
 
1204
         VTN(1) = IVX(K)
 
1205
         VTN(2) = IVY(K)
 
1206
         VPI    = VP(1)*IVX(K)+VP(TWO)*IVY(K)+VP(THREE)
 
1207
         VFU = 0.
 
1208
         DO 31 L = 1,NC
 
1209
            N    = L*3+1
 
1210
            N1   = L*4
 
1211
            DIX  = IVX(K)-VP(N1+1)
 
1212
            DIY  = IVY(K)-VP(N1+2)
 
1213
            EPES = DIX*DIX+DIY*DIY
 
1214
 
 
1215
            IF (BETA.LE.0) THEN
 
1216
               EP = EXP(EPES*FK(L))
 
1217
               CO = -VP(N1)*EP*2.*FK(L)
 
1218
            ELSE
 
1219
               EPE = 1.+EPES*FK(L)
 
1220
               EP  = EPE**(-BETA)
 
1221
               EP2 = EPE**(-BETA-1.)
 
1222
               CO  = VP(N1)*BETA*EP2*2*FK(L)
 
1223
            END IF
 
1224
 
 
1225
            VTN(N)   = EP
 
1226
            VTN(N+1) = CO*DIX
 
1227
            VTN(N+2) = CO*DIY
 
1228
            VFU      = VP(N1)*EP+VFU
 
1229
   31    CONTINUE
 
1230
 
 
1231
         VFU = VFU+VPI
 
1232
         DO 32 I = 1,NIN
 
1233
            VTP    = VTN(I)*WEI(K)
 
1234
            VC(I)  = VC(I)+(VZ(K)-VFU)*VTP
 
1235
            DO 36 J = 1,I
 
1236
               RMT(I,J) = RMT(I,J)+VTN(J)*VTP
 
1237
   36       CONTINUE
 
1238
   32    CONTINUE
 
1239
   30 CONTINUE
 
1240
 
 
1241
      DO 40 I = 2,NIN
 
1242
         DO 41 J = 1,I-1
 
1243
            RMT(J,I) = RMT(I,J)
 
1244
   41    CONTINUE
 
1245
   40 CONTINUE
 
1246
 
 
1247
      DO 50 I=1,NIN
 
1248
         RMT(I,I) = RMT(I,I)*(1+FS**2)
 
1249
   50 CONTINUE
 
1250
 
 
1251
      NCD = NIND
 
1252
      CALL LISIB(RMT,VC,NIN,NCD,SIG)
 
1253
 
 
1254
      IF (NCD.GT.0) THEN
 
1255
         DO 59 I = 1,3
 
1256
            VP(I) = VC(I)*VPES(I)+VP(I) 
 
1257
   59    CONTINUE
 
1258
         DO 60 I = 1,NC
 
1259
            DO 61 J = 4,6
 
1260
               INK     = J+4*(I-1)
 
1261
               INK1    = J+3*(I-1)
 
1262
               VP(INK) = VC(INK1)*VPES(J)+VP(INK)
 
1263
               IF(ABS(VP(INK)).GT.1000..AND.J.NE.4) NCD=-1
 
1264
   61        CONTINUE
 
1265
   60    CONTINUE
 
1266
 
 
1267
         IF (NCD.GT.0) THEN
 
1268
            SQM = 0.
 
1269
            DO 63 I = 1,NP
 
1270
               VG = VP(1)*IVX(I)+VP(TWO)*IVY(I)+VP(THREE)
 
1271
               DO 64 J = 1,NC
 
1272
                  N   = J*4
 
1273
                  EEE = ((VP(N+1)-IVX(I))**2 + (VP(N+2)-IVY(I))**2)/
 
1274
     2                                         (VP(N+3)**2)
 
1275
                  IF(BETA.LE.0.) THEN
 
1276
                     VG = VG+VP(N)*EXP(-EEE*4*ALOG(2.))
 
1277
                  ELSE
 
1278
                     VG = VG+VP(N)*(1+EEE)**(-BETA)
 
1279
                  END IF
 
1280
   64          CONTINUE
 
1281
               SQM = SQM+(VZ(I)-VG)**2*WEI(I)
 
1282
   63       CONTINUE
 
1283
            SQM = SQM/(NP-NIN)
 
1284
         END IF
 
1285
      END IF
 
1286
 
 
1287
      IF(NCD.LT.1) LG=1
 
1288
      RETURN
 
1289
      END