← Back to branch summary

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

« back to all changes in this revision

Viewing changes to stdred/spec/libsrc/prc.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/spec/libsrc/prc.for
 
1
C @(#)prc.for   19.1 (ES0-DMD) 02/25/03 14:29:02
 
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
CC
 
30
CC  *STATIS* : STANDARDIZATION OF THE OBSERVATIONS AND
 
31
CC               CALCULATION OF SOME STATISTICS
 
32
CC
 
33
       SUBROUTINE STATIS(X,Y,XMED,XMAD,RESDU,WEIGHTS,AW,
 
34
     1 NMVAL,C,H,JNDVC,NSTOP,NVAR,NVAD,NVSB,NCAS,JCST,JPRT,
 
35
     1 NDXX,NDXY,LUB,NMXV,PREC,LAB)
 
36
 
 
37
       INCLUDE 'MID_REL_INCL:implicit.inc'
 
38
 
 
39
       INTEGER J,JCST,JNC,JPRT,JS,JSM,JSPA,JSPB,JSPP,LUB
 
40
       INTEGER NCAS,NDXX,NDXY,NMXV,NVAD,NSTOP,NVAR,NVSB
 
41
       REAL AL,AMDAN,DIAG,SPEAR,PEARC,PREC
 
42
C       
 
43
       REAL X(NDXX,NDXY),Y(NDXY),XMED(NDXX),XMAD(NDXX)
 
44
       REAL RESDU(NDXY),WEIGHTS(NDXY),AW(NDXY),C(NMXV,NDXX)
 
45
       INTEGER NMVAL(NDXY),JNDVC(NDXX)       
 
46
       
 
47
C       DIMENSION X(NDXX,NDXY),Y(NDXY),XMED(NDXX),XMAD(NDXX)
 
48
C       DIMENSION RESDU(NDXY),WEIGHTS(NDXY),AW(NDXY),NMVAL(NDXY)
 
49
C       DIMENSION C(NMXV,NDXX),JNDVC(NDXX)
 
50
 
 
51
       DOUBLE PRECISION H(NMXV,NDXX)
 
52
       CHARACTER*10 LAB(NDXX)
 
53
CC-----IN THIS SUBROUTINE THE VECTOR "RESDU" WILL BE USED AS A
 
54
CC-----WORKING VECTOR
 
55
       AL=NCAS
 
56
       IF (JCST.NE.0) GOTO 60
 
57
       DO 50 J=1,NVAD
 
58
       XMED(J)=0.0
 
59
       DO 10 JNC=1,NCAS
 
60
 10    RESDU(JNC)=ABS(X(J,JNC))
 
61
       XMAD(J)=AMDAN(AW,NDXY,RESDU,NCAS)*1.4826
 
62
       IF (ABS(XMAD(J)).GT.PREC) GOTO 30
 
63
       XMAD(J)=0.0
 
64
       DO 20 JNC=1,NCAS
 
65
 20    XMAD(J)=XMAD(J)+RESDU(JNC)
 
66
       XMAD(J)=(XMAD(J)/AL)*1.2533
 
67
       IF (ABS(XMAD(J)).GT.PREC) GOTO 30
 
68
       WRITE(LUB,8000) LAB(J)
 
69
       NSTOP=1
 
70
       RETURN
 
71
 30    DO 40 JNC=1,NCAS
 
72
       X(J,JNC)=X(J,JNC)/XMAD(J)
 
73
 40    CONTINUE
 
74
 50    CONTINUE
 
75
       WRITE(LUB,8015)
 
76
       GOTO 150
 
77
 60    DO 120 J=1,NVAD
 
78
       IF(J.EQ.NVAR) GOTO 120
 
79
       DO 70 JNC=1,NCAS
 
80
 70    RESDU(JNC)=X(J,JNC)
 
81
       XMED(J)=AMDAN(AW,NDXY,RESDU,NCAS)
 
82
       DO 80 JNC=1,NCAS
 
83
 80    RESDU(JNC)=ABS(RESDU(JNC)-XMED(J))
 
84
       XMAD(J)=AMDAN(AW,NDXY,RESDU,NCAS)*1.4826
 
85
       IF (ABS(XMAD(J)).GT.PREC) GOTO 100
 
86
       XMAD(J)=0.0
 
87
       DO 90 JNC=1,NCAS
 
88
 90    XMAD(J)=XMAD(J)+RESDU(JNC)
 
89
       XMAD(J)=(XMAD(J)/AL)*1.2533
 
90
       IF (ABS(XMAD(J)).GT.PREC) GOTO 100
 
91
       WRITE(LUB,8000) LAB(J)
 
92
       NSTOP=1
 
93
       RETURN
 
94
 100   DO 110 JNC=1,NCAS
 
95
       X(J,JNC)=(X(J,JNC)-XMED(J))/XMAD(J)
 
96
 110   CONTINUE
 
97
 120   CONTINUE
 
98
       XMED(NVAR)=1.0
 
99
       WRITE(LUB,8010)
 
100
       IF (JCST.EQ.1) WRITE(LUB,8020) (LAB(J),J=1,NVSB),LAB(NVAD)
 
101
       IF (JCST.EQ.0) WRITE(LUB,8020) (LAB(J),J=1,NVAD)
 
102
       IF ((NVAD.GT.6.AND.JCST.EQ.0).OR.(NVAD.GT.7.AND.JCST.EQ.1))
 
103
     1 GOTO 130
 
104
       IF (JCST.EQ.0) WRITE(LUB,8030)(XMED(J),J=1,NVAD)
 
105
       IF (JCST.EQ.1) WRITE(LUB,8030)(XMED(J),J=1,NVSB),XMED(NVAD)
 
106
       GOTO 140
 
107
 130   WRITE(LUB,8030)(XMED(J),J=1,6)
 
108
       IF (JCST.EQ.0) WRITE(LUB,8040)(XMED(J),J=7,NVAD)
 
109
       IF (JCST.EQ.1.AND.NVSB.GE.7)
 
110
     1 WRITE(LUB,8040)(XMED(J),J=7,NVSB),XMED(NVAD)
 
111
       IF (JCST.EQ.1.AND.NVSB.LT.7) WRITE(LUB,8040) XMED(NVAD)
 
112
 140   XMED(NVAR)=0.0
 
113
       WRITE(LUB,8050)
 
114
       IF (JCST.EQ.1) WRITE(LUB,8020) (LAB(J),J=1,NVSB),LAB(NVAD)
 
115
       IF (JCST.EQ.0) WRITE(LUB,8020) (LAB(J),J=1,NVAD)
 
116
 150   IF((NVAD.GT.6.AND.JCST.EQ.0).OR.(NVAD.GT.7.AND.JCST.EQ.1))
 
117
     1 GOTO 160
 
118
       IF (JCST.EQ.0) WRITE(LUB,8030)(XMAD(J),J=1,NVAD)
 
119
       IF (JCST.EQ.1) WRITE(LUB,8030)(XMAD(J),J=1,NVSB),XMAD(NVAD)
 
120
       GOTO 170
 
121
 160   WRITE(LUB,8030)(XMAD(J),J=1,6)
 
122
       IF (JCST.EQ.0) WRITE(LUB,8040)(XMAD(J),J=7,NVAD)
 
123
       IF (JCST.EQ.1.AND.NVSB.GE.7)
 
124
     1 WRITE(LUB,8040)(XMAD(J),J=7,NVSB),XMAD(NVAD)
 
125
       IF (JCST.EQ.1.AND.NVSB.LT.7) WRITE(LUB,8040) XMAD(NVAD)
 
126
 170   IF(JCST.NE.0) XMAD(NVAR)=1.0
 
127
       IF (JPRT.LT.2) GOTO 200
 
128
       WRITE(LUB,8060)
 
129
       IF (JCST.EQ.1) WRITE(LUB,8020) (LAB(J),J=1,NVSB),LAB(NVAD)
 
130
       IF (JCST.EQ.0) WRITE(LUB,8020) (LAB(J),J=1,NVAD)
 
131
       DO 190 JNC=1,NCAS
 
132
       IF((NVAD.GT.6.AND.JCST.EQ.0).OR.(NVAD.GT.7.AND.JCST.EQ.1))
 
133
     1 GOTO 180
 
134
       IF (JCST.EQ.1) WRITE(LUB,8070) NMVAL(JNC),(X(J,JNC),J=1,NVSB),
 
135
     1 X(NVAD,JNC)
 
136
       IF (JCST.EQ.0) WRITE(LUB,8070) NMVAL(JNC),(X(J,JNC),J=1,NVAD)
 
137
       GOTO 190
 
138
 180   WRITE(LUB,8070) NMVAL(JNC),(X(J,JNC),J=1,6)
 
139
       IF (JCST.EQ.1.AND.NVSB.GE.7)
 
140
     1 WRITE(LUB,8040) (X(J,JNC),J=7,NVSB),X(NVAD,JNC)
 
141
       IF (JCST.EQ.1.AND.NVSB.LT.7) WRITE(LUB,8040) X(NVAD,JNC)
 
142
       IF (JCST.EQ.0) WRITE(LUB,8040) (X(J,JNC),J=7,NVAD)
 
143
 190   CONTINUE
 
144
 200   IF(JPRT.EQ.0) GOTO 260
 
145
       DO 220 JSPA=1,NVAR
 
146
       IF(JSPA.EQ.NVAR.AND.JCST.EQ.1) GOTO 220
 
147
       JS=JSPA+1
 
148
       DO 210 JSPB=JS,NVAD
 
149
       IF(JSPB.EQ.NVAR.AND.JCST.EQ.1) GOTO 210
 
150
       JSM=JSPB-1
 
151
       CALL CORR(X,JSPB,JSPA,NCAS,NVSB,JCST,AW,Y,RESDU,NDXX,NDXY,
 
152
     1 SPEAR,PEARC,PREC)
 
153
       C(JSM,JSPA)=SPEAR
 
154
       H(JSM,JSPA)=PEARC
 
155
 210   CONTINUE
 
156
 220   CONTINUE
 
157
       WRITE(LUB,8080) LAB(NVAD)
 
158
       DIAG=1.0
 
159
       WRITE(LUB,8090) LAB(1),DIAG
 
160
       DO 230 JSPA=1,NVSB
 
161
       JSPP=JSPA+1
 
162
       IF(JCST.EQ.1.AND.JSPP.EQ.NVAR) GOTO 230
 
163
       WRITE(LUB,8090) LAB(JSPP),(H(JSPA,JSPB),JSPB=1,JSPA),DIAG
 
164
 230   CONTINUE
 
165
       IF(JCST.EQ.1) WRITE(LUB,8090) LAB(NVAD),(H(JSPA,JSPB),
 
166
     1 JSPB=1,NVSB),DIAG
 
167
       IF (JCST.EQ.0) WRITE(LUB,8090) LAB(NVAD),(H(JSPA,JSPB),
 
168
     1 JSPB=1,NVAR),DIAG
 
169
       DO 240 J=1,NVAR
 
170
 240   JNDVC(J)=J
 
171
       WRITE(LUB,8100) LAB(NVAD)
 
172
       WRITE(LUB,8090) LAB(1),DIAG
 
173
       DO 250 JSPA=1,NVSB
 
174
       JSPP=JSPA+1
 
175
       IF(JCST.EQ.1.AND.JSPP.EQ.NVAR) GOTO 250
 
176
       WRITE(LUB,8090) LAB(JSPP),(C(JSPA,JSPB),JSPB=1,JSPA),DIAG
 
177
 250   CONTINUE
 
178
       IF(JCST.EQ.1) WRITE(LUB,8090) LAB(NVAD),(C(JSPA,JSPB),
 
179
     1 JSPB=1,NVSB),DIAG
 
180
       IF (JCST.EQ.0) WRITE(LUB,8090) LAB(NVAD),(C(JSPA,JSPB),
 
181
     1 JSPB=1,NVAR),DIAG
 
182
 260   DO 270 JNC=1,NCAS
 
183
CC-----MEANWHILE, INITIALIZATION OF THE WEIGHTS
 
184
       WEIGHTS(JNC)=1.0
 
185
       Y(JNC)=X(NVAD,JNC)
 
186
 270   CONTINUE
 
187
 8000  FORMAT(/' VARIABLE WITH LABEL ',A10,' IS CONSTANT OVER'/
 
188
     1 ' ALL OBSERVATIONS HENCE THE PROGRAM CANNOT RUN.',
 
189
     1 /' PLEASE REPEAT THE ANALYSIS WITHOUT THIS VARIABLE.'/)
 
190
 8010  FORMAT(//' MEDIANS  = '/)
 
191
 8015  FORMAT(//' DISPERSION OF ABSOLUTE VALUES = '/)
 
192
 8020  FORMAT(7X,6(A10,1X)/15X,5(A10,1X))
 
193
 8030  FORMAT(7X,6(F10.4,1X))
 
194
 8040  FORMAT(15X,5(F10.4,1X))
 
195
 8050  FORMAT(//' DISPERSIONS = '/)
 
196
 8060  FORMAT(//' THE STANDARDIZED OBSERVATIONS ARE:'/)
 
197
 8070  FORMAT(I5,2X,6(F10.4,1X))
 
198
 8080  FORMAT(//' PEARSON CORRELATION COEFFICIENTS',
 
199
     1 ' BETWEEN THE VARIABLES'/1X,'( ',
 
200
     1 A10,' IS THE RESPONSE VARIABLE)'/)
 
201
 8090  FORMAT(1X,A10,(2X,11F6.2)/)
 
202
 8100  FORMAT(//' SPEARMAN RANK CORRELATION COEFFICIENTS',
 
203
     1 ' BETWEEN THE VARIABLES'/1X,'( ',A10,
 
204
     1 ' IS THE RESPONSE VARIABLE)'/ )
 
205
       RETURN
 
206
       END
 
207
CC
 
208
CC  *CORR* : CALCULATES THE SPEARMAN RANK CORRELATION COEFFICIENT
 
209
CC             BETWEEN VARIABLES JSPA AND JSPB.
 
210
       SUBROUTINE CORR(X,JSPB,JSPA,NCAS,NVSB,JCST,AW,Y,RESDU,NDXX,NDXY,
 
211
     1 SPEAR,PEARC,PREC)
 
212
 
 
213
C       INCLUDE 'MID_REL_INCL:implicit.inc'
 
214
C       DIMENSION X(NDXX,NDXY),AW(NDXY),Y(NDXY),RESDU(NDXY)
 
215
        INTEGER JCST,JSPA,JSPB,JTA,K,KK,NCM,NDX,NDXX,NDXY,NCAS
 
216
        INTEGER NJ,NJA,NJAP,NJB,NJP,NTWE,NVSB
 
217
        REAL AL,CVVAR,PEARC,PREC,SPEAR,SUMA,SUMB,STDA,STDB,TA,TB,W
 
218
 
 
219
        REAL X(NDXX,NDXY),AW(NDXY),Y(NDXY),RESDU(NDXY)
 
220
       
 
221
       
 
222
CC-----IN THIS SUBROUTINE THE VECTOR "RESDU" IS USED AS A WORKING
 
223
CC-----VECTOR FOR CALCULATING THE CORRELATION MATRICES
 
224
       SPEAR=0.0
 
225
       SUMA=0.0
 
226
       SUMB=0.0
 
227
       STDA=0.0
 
228
       STDB=0.0
 
229
       CVVAR=0.0
 
230
       AL=NCAS
 
231
       IF (((JSPB-JSPA).EQ.1).OR.(JCST.EQ.1.AND.JSPA.EQ.NVSB))
 
232
     1 GOTO  5
 
233
       GOTO 105
 
234
 5     NTWE=0
 
235
       DO 10 NJ=1,NCAS
 
236
       AW(NJ)=X(JSPA,NJ)
 
237
       RESDU(NJ)=NJ
 
238
 10    CONTINUE
 
239
       NCM=NCAS-1
 
240
 15    DO 20 NJ=1,NCM
 
241
       NJP=NJ+1
 
242
       NJA=NJ
 
243
       DO 30 NJB=NJP,NCAS
 
244
       IF(AW(NJB).GE.AW(NJA)) GOTO 30
 
245
       NJA=NJB
 
246
 30    CONTINUE
 
247
       IF(NJA.EQ.NJ) GOTO 20
 
248
       W=AW(NJ)
 
249
       AW(NJ)=AW(NJA)
 
250
       AW(NJA)=W
 
251
       W=RESDU(NJA)
 
252
       RESDU(NJA)=RESDU(NJ)
 
253
       RESDU(NJ)=W
 
254
 20    CONTINUE
 
255
       DO 40 NJ=1,NCAS
 
256
       NDX=RESDU(NJ)
 
257
       IF(NTWE.EQ.0) Y(NDX)=NJ
 
258
       IF(NTWE.NE.0) AW(NDX)=NJ
 
259
 40    CONTINUE
 
260
       NJ=1
 
261
 55    TA=1
 
262
       NJA=RESDU(NJ)
 
263
       IF(NTWE.EQ.0) TB=Y(NJA)
 
264
       IF(NTWE.NE.0) TB=AW(NJA)
 
265
 60    NJA=RESDU(NJ)
 
266
       NJAP=RESDU(NJ+1)
 
267
       IF(NTWE.EQ.0.AND.(X(JSPA,NJA).NE.X(JSPA,NJAP))) GOTO 70
 
268
       IF(NTWE.NE.0.AND.(X(JSPB,NJA).NE.X(JSPB,NJAP))) GOTO 70
 
269
       TA=TA+1.0
 
270
       IF(NTWE.EQ.0) TB=TB+Y(NJAP)
 
271
       IF(NTWE.NE.0) TB=TB+AW(NJAP)
 
272
       NJ=NJ+1
 
273
       IF (NJ.EQ.NCAS) GOTO 70
 
274
       GOTO 60
 
275
 70    IF(TA.EQ.1.0) GOTO 100
 
276
       JTA=TA
 
277
       NJB=NJ+1-JTA
 
278
       DO 80 K=NJB,NJ
 
279
       KK=RESDU(K)
 
280
       IF(NTWE.EQ.0) Y(KK)=TB/TA
 
281
       IF(NTWE.NE.0) AW(KK)=TB/TA
 
282
 80    CONTINUE
 
283
 100   NJ=NJ+1
 
284
       IF (NJ.LT.NCAS) GOTO 55
 
285
       IF(NTWE.NE.0) GOTO 150
 
286
 105   DO 110 NJ=1,NCAS
 
287
       AW(NJ)=X(JSPB,NJ)
 
288
       RESDU(NJ)=NJ
 
289
       SUMA=SUMA+X(JSPA,NJ)
 
290
       SUMB=SUMB+X(JSPB,NJ)
 
291
 110   CONTINUE
 
292
       SUMA=SUMA/AL
 
293
       SUMB=SUMB/AL
 
294
       NTWE=1
 
295
       GOTO 15
 
296
 150   DO 120 NJ=1,NCAS
 
297
       SPEAR=SPEAR+(Y(NJ)-AW(NJ))**2
 
298
       STDA=STDA+(X(JSPA,NJ)-SUMA)**2
 
299
       STDB=STDB+(X(JSPB,NJ)-SUMB)**2
 
300
       CVVAR=CVVAR+(X(JSPA,NJ)-SUMA)*(X(JSPB,NJ)-SUMB)
 
301
 120   CONTINUE
 
302
       SPEAR=1.0-6.0*SPEAR/(AL*(AL**2-1.0))
 
303
       STDA=SQRT(STDA*STDB)
 
304
       IF(ABS(STDA).GT.PREC) PEARC=CVVAR/(STDA)
 
305
       IF(ABS(STDA).LE.PREC) PEARC=99.99
 
306
       RETURN
 
307
       END
 
308
CC
 
309
CC  *PULL* : SEARCHES THE KTH VALUE IN A VECTOR
 
310
CC             OF LENGTH 'NCAS'.
 
311
       FUNCTION PULL(AW,NDXY,AA,NCAS,K)
 
312
 
 
313
C       INCLUDE 'MID_REL_INCL:implicit.inc'
 
314
C       DIMENSION AA(NCAS),AW(NDXY)
 
315
 
 
316
       INTEGER J,JNC,K,L,LR,NCAS,NDXY
 
317
       REAL AX,PULL,WA
 
318
 
 
319
       REAL AA(NCAS),AW(NDXY)
 
320
 
 
321
 
 
322
 
 
323
 
 
324
       DO 10 JNC=1,NCAS
 
325
 10    AW(JNC)=AA(JNC)
 
326
       L=1
 
327
       LR=NCAS
 
328
 20    IF(L.GE.LR) GOTO 90
 
329
       AX=AW(K)
 
330
       JNC=L
 
331
       J=LR
 
332
 30    IF(JNC.GT.J) GOTO 80
 
333
 40    IF(AW(JNC).GE.AX) GOTO 50
 
334
       JNC=JNC+1
 
335
       GOTO 40
 
336
 50    IF(AW(J).LE.AX) GOTO 60
 
337
       J=J-1
 
338
       GOTO 50
 
339
 60    IF(JNC.GT.J) GOTO 70
 
340
       WA=AW(JNC)
 
341
       AW(JNC)=AW(J)
 
342
       AW(J)=WA
 
343
       JNC=JNC+1
 
344
       J=J-1
 
345
 70    GOTO 30
 
346
 80    IF(J.LT.K) L=JNC
 
347
       IF (K.LT.JNC) LR=J
 
348
       GOTO 20
 
349
 90    PULL=AW(K)
 
350
       RETURN
 
351
       END
 
352
CC
 
353
CC  *AMDAN* : CALCULATES THE MEDIAN OF A VECTOR
 
354
CC              OF LENGTH 'NCAS'.
 
355
 
 
356
       FUNCTION AMDAN(AW,NDXY,AA,NCAS)
 
357
 
 
358
C       INCLUDE 'MID_REL_INCL:implicit.inc'
 
359
C       DIMENSION AA(NCAS),AW(NDXY)
 
360
 
 
361
       INTEGER JNDL,NCAS,NDXY
 
362
       REAL AL,AMDAN,PULL
 
363
 
 
364
       REAL AA(NCAS),AW(NDXY)
 
365
 
 
366
       AL=NCAS
 
367
       JNDL=INT(AL/2.0)
 
368
       IF (MOD(NCAS,2).EQ.0) THEN
 
369
          AMDAN= (PULL(AW,NDXY,AA,NCAS,JNDL)+
 
370
     1                    PULL(AW,NDXY,AA,NCAS,JNDL+1))/2.0
 
371
                             ELSE
 
372
       AMDAN= PULL(AW,NDXY,AA,NCAS,JNDL+1)
 
373
       ENDIF
 
374
       RETURN
 
375
       END
 
376
CC
 
377
CC  *RSQU* : CALCULATES COEFFICIENT OF DETERMINATION AND F-VALUE
 
378
CC           OF LS OR RLS
 
379
CC
 
380
       FUNCTION RSQU(NCAS,NVAD,JCST,Y,NDXY,SSE,FVALUE,PREC,
 
381
     1 XMAD,XMED,NDXX,WEIGHTS,NNUL)
 
382
 
 
383
C       INCLUDE 'MID_REL_INCL:implicit.inc'
 
384
C       DIMENSION Y(NDXY),XMAD(NDXX),XMED(NDXX),WEIGHTS(NDXY)
 
385
 
 
386
       INTEGER  JCST,JNC,NCAS,NDXX,NDXY,NNUL,NVAD
 
387
       REAL  AL,DF1,DF2,FVALUE,PREC,RSQU,RYM,SSE,SSEL,SST
 
388
 
 
389
       REAL Y(NDXY),XMAD(NDXX),XMED(NDXX),WEIGHTS(NDXY)
 
390
 
 
391
 
 
392
 
 
393
       RYM=0.0
 
394
       IF (JCST.EQ.0) GOTO 20
 
395
       DO 10 JNC=1,NCAS
 
396
 10    RYM=(Y(JNC)*XMAD(NVAD)+XMED(NVAD))*WEIGHTS(JNC)+RYM
 
397
       AL=NNUL
 
398
       RYM=RYM/AL
 
399
 20    SST=0.0
 
400
       DO 30 JNC=1,NCAS
 
401
 30    SST=SST+(((Y(JNC)*XMAD(NVAD)+XMED(NVAD))-RYM)**2)*WEIGHTS(JNC)
 
402
       DF1=(NVAD-1)-JCST
 
403
       DF2=NNUL-(NVAD-1)
 
404
       IF(SST.LT.PREC)SST=PREC
 
405
       RSQU=1.0-(SSE/SST)
 
406
       IF (RSQU.LT.0.0) RSQU=0.0
 
407
       IF (RSQU.GT.1.0) RSQU=1.0
 
408
       SSEL=SSE
 
409
       IF(SSEL.LT.PREC)SSEL=PREC
 
410
       FVALUE=((SST-SSEL)/DF1)/(SSEL/DF2)
 
411
       IF(FVALUE.LT.(0.0))FVALUE=0.0
 
412
       RETURN
 
413
       END
 
414
CC
 
415
CC  *RESTT* : CALCULATES FOR ALL JNC (JNC=1,..,NCAS) THE RESIDUAL
 
416
CC              Y(JNC)-SUM (GESCF(J)*X(J,JNC)), FOR J=1,..,LTSTE.
 
417
 
 
418
       SUBROUTINE RESTT(GESCF,JABS,JTR,LTSTE,NCAS,NVAD,NZWE,
 
419
     1 X,Y,RESDU,WEIGHTS,XMED,XMAD,NDXX,NDXY,NOPT,SCAL,QUAN,PREC)
 
420
 
 
421
C       INCLUDE 'MID_REL_INCL:implicit.inc'
 
422
C       DIMENSION GESCF(LTSTE),X(NDXX,NDXY),Y(NDXY),WEIGHTS(NDXY)
 
423
C       DIMENSION RESDU(NDXY),XMED(NDXX),XMAD(NDXX)
 
424
 
 
425
       INTEGER J,JABS,JNC,JTR,LTSTE,NCAS,NDXX,NDXY,NOPT,NVAD,NZWE
 
426
       REAL ANVAR,AVLQS,PREC,QUAN,SCAL,WSC
 
427
 
 
428
       REAL GESCF(LTSTE),X(NDXX,NDXY),Y(NDXY),WEIGHTS(NDXY)
 
429
       REAL RESDU(NDXY),XMED(NDXX),XMAD(NDXX)
 
430
 
 
431
 
 
432
 
 
433
CC-----JTR = 1, THEN CALCULATE THE RESIDUALS FOR THE STANDARDIZED OBSERVATIONS
 
434
CC-----      0         "      "             "              RAW
 
435
CC-----JABS = 1, CALCULATE THE ABSOLUTE RESIDUALS
 
436
CC-----       0                RESIDUALS
 
437
CC-----NOPT = 1, CALCULATE THE LMS FINAL SCALE AND THE RESULTING WEIGHTS
 
438
CC-----       0, OTHERWISE
 
439
       NZWE=0
 
440
       SCAL=0.0
 
441
       AVLQS=0.0
 
442
       ANVAR=NVAD-1
 
443
       DO 10 JNC=1,NCAS
 
444
       IF (JTR.NE.1)  THEN
 
445
          RESDU(JNC)=Y(JNC)
 
446
                      ELSE
 
447
          RESDU(JNC)=Y(JNC)*XMAD(NVAD)+XMED(NVAD)
 
448
       ENDIF
 
449
       DO 20 J=1,LTSTE
 
450
       IF (JTR.NE.1) THEN
 
451
          RESDU(JNC)=RESDU(JNC)-X(J,JNC)*GESCF(J)
 
452
                     ELSE
 
453
          RESDU(JNC)=RESDU(JNC)-((X(J,JNC)*XMAD(J)+XMED(J))*GESCF(J))
 
454
       ENDIF
 
455
 20    CONTINUE
 
456
       IF (JABS.EQ.1.AND.NOPT.NE.1) RESDU(JNC)=ABS(RESDU(JNC))
 
457
       IF (ABS(RESDU(JNC)).LT.PREC) NZWE=NZWE+1
 
458
       IF (NOPT.EQ.0) GOTO 10
 
459
       IF (ABS(RESDU(JNC)).LE.(2.5*QUAN)) THEN
 
460
          WSC=1.0
 
461
                                     ELSE
 
462
          WSC=0.0
 
463
       ENDIF
 
464
       AVLQS=AVLQS+WSC
 
465
       SCAL=SCAL+(RESDU(JNC)*RESDU(JNC))*WSC
 
466
 10    CONTINUE
 
467
C       write(6,*) 'In PRC/RESTT: NOPT = ',NOPT
 
468
       IF (NOPT.EQ.0) RETURN
 
469
       SCAL=SQRT(SCAL/(AVLQS-ANVAR)) + PREC
 
470
       NZWE=0
 
471
       DO 40 JNC=1,NCAS
 
472
       IF (ABS(RESDU(JNC)).LE.(2.5*SCAL)) THEN
 
473
          WEIGHTS(JNC)=1.0
 
474
          NZWE=NZWE+1
 
475
                                          ELSE
 
476
          WEIGHTS(JNC)=0.0
 
477
       ENDIF
 
478
C       write(6,*) 'In PRC: ',jnc,resdu(jnc),scal,nzwe,weights(jnc)
 
479
 40    CONTINUE
 
480
       RETURN
 
481
       END