]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 1 | * |
2 | * $Id$ | |
3 | * | |
4 | * $Log$ | |
5 | * Revision 1.1.1.1 1996/04/01 15:02:28 mclareni | |
6 | * Mathlib gen | |
7 | * | |
8 | * | |
9 | #include "gen/pilot.h" | |
10 | #if defined(CERNLIB_CDC) | |
11 | SUBROUTINE PARLSQ | |
12 | 1 (PARX2P,PARY2P,NPAR2P,COEF2P,SDEV2P) | |
13 | C | |
14 | C TO FIT A PARABOLA TO NPAR2P POINTS | |
15 | C | |
16 | C NPAR2P NO. OF POINTS | |
17 | C PARX2P(I) X VALUE OF POINT I | |
18 | C PARY2P(I) Y VALUE OF POINT I | |
19 | C | |
20 | C COEF2P(1...3) COEFFICIENTS OF THE FITTED PARABOLA | |
21 | C Y=COEF2P(1) + COEF2P(2)*X + COEF2P(3)*X**2 | |
22 | C SDEV2P= VARIANCE | |
23 | C METHOD : CHI**2 = MIN EQUATION SOLVED EXPLICITLY | |
24 | DIMENSION PARX2P(NPAR2P),PARY2P(NPAR2P),COEF2P(NPAR2P) | |
25 | C | |
26 | DO 3 I=1,3 | |
27 | 3 COEF2P(I)=0. | |
28 | SDEV2P=0. | |
29 | IF(NPAR2P.LT.3) GO TO 10 | |
30 | F=NPAR2P | |
31 | C--- CENTER X VALUES FOR REASONS OF MACHINE PRECISION | |
32 | XM=0. | |
33 | DO 2 I=1,NPAR2P | |
34 | 2 XM=XM+PARX2P(I) | |
35 | XM=XM/F | |
36 | X2=0. | |
37 | X3=0. | |
38 | X4=0. | |
39 | Y=0. | |
40 | Y2=0. | |
41 | XY=0. | |
42 | X2Y=0. | |
43 | DO 1 I=1,NPAR2P | |
44 | S=PARX2P(I)-XM | |
45 | T=PARY2P(I) | |
46 | S2=S*S | |
47 | X2=X2+S2 | |
48 | X3=X3+S*S2 | |
49 | X4=X4+S2*S2 | |
50 | Y=Y+T | |
51 | Y2=Y2+T*T | |
52 | XY=XY+S*T | |
53 | X2Y=X2Y+S2*T | |
54 | 1 CONTINUE | |
55 | A=(F*X4-X2**2)*X2-F*X3**2 | |
56 | IF(A.EQ.0.) GOTO 10 | |
57 | COEF2P(3)=(X2*(F*X2Y-X2*Y)-F*X3*XY)/A | |
58 | COEF2P(2)=(XY-X3*COEF2P(3))/X2 | |
59 | COEF2P(1)=(Y-X2*COEF2P(3))/F | |
60 | IF(NPAR2P.EQ.3) GOTO 6 | |
61 | SDEV2P=Y2-(COEF2P(1)*Y+COEF2P(2)*XY+COEF2P(3)*X2Y) | |
62 | IF(SDEV2P.LT.0.) SDEV2P=0. | |
63 | SDEV2P=SDEV2P/(F-3) | |
64 | 6 COEF2P(1)=COEF2P(1)+XM*(XM*COEF2P(3)-COEF2P(2)) | |
65 | COEF2P(2)=COEF2P(2)-2.*XM*COEF2P(3) | |
66 | 10 CONTINUE | |
67 | RETURN | |
68 | END | |
69 | #endif | |
70 | #if !defined(CERNLIB_CDC) | |
71 | SUBROUTINE PARLSQ | |
72 | 1 (PARX2P,PARY2P,NPAR2P,COEF2P,SDEV2P) | |
73 | C | |
74 | C TO FIT A PARABOLA TO NPAR2P POINTS | |
75 | C | |
76 | C NPAR2P NO. OF POINTS | |
77 | C PARX2P(I) X VALUE OF POINT I | |
78 | C PARY2P(I) Y VALUE OF POINT I | |
79 | C | |
80 | C COEF2P(1...3) COEFFICIENTS OF THE FITTED PARABOLA | |
81 | C Y=COEF2P(1) + COEF2P(2)*X + COEF2P(3)*X**2 | |
82 | C SDEV2P= VARIANCE | |
83 | C METHOD : CHI**2 = MIN EQUATION SOLVED EXPLICITLY | |
84 | DIMENSION PARX2P(NPAR2P),PARY2P(NPAR2P),COEF2P(NPAR2P) | |
85 | DOUBLE PRECISION A,CZ(3),F,S,S2,T,XM,XY,X2,X2Y,X3,X4,Y,Y2 | |
86 | C | |
87 | DO 3 I=1,3 | |
88 | 3 CZ(I)=0.D0 | |
89 | SDEV2P=0. | |
90 | IF(NPAR2P.LT.3) GO TO 10 | |
91 | F=NPAR2P | |
92 | C--- CENTER X VALUES FOR REASONS OF MACHINE PRECISION | |
93 | XM=0.D0 | |
94 | DO 2 I=1,NPAR2P | |
95 | 2 XM=XM+PARX2P(I) | |
96 | XM=XM/F | |
97 | X2=0.D0 | |
98 | X3=0.D0 | |
99 | X4=0.D0 | |
100 | Y=0.D0 | |
101 | Y2=0.D0 | |
102 | XY=0.D0 | |
103 | X2Y=0.D0 | |
104 | DO 1 I=1,NPAR2P | |
105 | S=PARX2P(I)-XM | |
106 | T=PARY2P(I) | |
107 | S2=S*S | |
108 | X2=X2+S2 | |
109 | X3=X3+S*S2 | |
110 | X4=X4+S2*S2 | |
111 | Y=Y+T | |
112 | Y2=Y2+T*T | |
113 | XY=XY+S*T | |
114 | X2Y=X2Y+S2*T | |
115 | 1 CONTINUE | |
116 | A=(F*X4-X2**2)*X2-F*X3**2 | |
117 | IF(A.EQ.0.D0) GOTO 10 | |
118 | CZ(3)=(X2*(F*X2Y-X2*Y)-F*X3*XY)/A | |
119 | CZ(2)=(XY-X3*CZ(3))/X2 | |
120 | CZ(1)=(Y-X2*CZ(3))/F | |
121 | IF(NPAR2P.EQ.3) GOTO 6 | |
122 | SDEV2P=Y2-(CZ(1)*Y+CZ(2)*XY+CZ(3)*X2Y) | |
123 | IF(SDEV2P.LT.0.) SDEV2P=0. | |
124 | SDEV2P=SDEV2P/(F-3.D0) | |
125 | 6 CZ(1)=CZ(1)+XM*(XM*CZ(3)-CZ(2)) | |
126 | CZ(2)=CZ(2)-2.*XM*CZ(3) | |
127 | 10 CONTINUE | |
128 | DO 11 I=1,3 | |
129 | 11 COEF2P(I)=CZ(I) | |
130 | RETURN | |
131 | END | |
132 | #endif |