[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / e / e406cod.inc

*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1996/04/01 15:02:28  mclareni
* Mathlib gen
*
*
*
* e406cod.inc
*
      DIMENSION C(0:128),CC(0:128),AC(0:64),BC(0:64)

      PARAMETER (Z1 = 1, Z2 = 2, HF = Z1/2)

      ALFA=HF*(B-A)
      BETA=HF*(B+A)
      C1=F(A)
      C2=F(B)
      AC(0)=C2+C1
      AC(1)=C2-C1

      DO 1 I = 1,7
      I1=2**(I-1)
      I2=I1-1
      I3=2*I1
      C1=Z2/I1
      C2=PI/I1

      DO 2 J = 0,I2
    2 C(J)=F(ALFA*COS((J+HF)*C2)+BETA)

C    COMPUTE B-COEFFICIENTS

      DO 3 J = 0,I2
      F1=J*C2
      F2=-HF*F1

      C3=2*COS(F1)
      A2=0
      A1=0
      A0=C(I2)
      DO 4 K = I2-1,0,-1
      A2=A1
      A1=A0
    4 A0=C(K)+C3*A1-A2
    3 BC(J)=C1*(A0*COS(F1+F2)-A1*COS(F2))
      BC(I1)=0

C    COMPUTE NEW C-COEFFICIENTS

      DO 5 J = 0,I1
      C(J)=HF*(AC(J)+BC(J))
    5 C(I1+J)=HF*(AC(I1-J)-BC(I1-J))

C    CHECK IF CALCULATION IS FINISHED

      CMX=0
      DO 6 J = 0,I3
      CC(J)=ABS(C(J))
    6 CMX=MAX(CMX,CC(J))

      IF(CMX .GT. 0) THEN
       CMX=1/CMX
       CC(I3)=HF*CC(I3)
       A0=CC(I2)*CMX
       A1=CC(I1)*CMX
       DO 10 J = I1+2,I3
       A2=CC(J)*CMX
       IF(A0 .LE. EPS .AND. A1. LE. EPS .AND. A2 .LE. EPS) GO TO 9

       A0=A1
   10  A1=A2
      ENDIF

C    DOUBLE THE NUMBER OF COEFFICIENTS.

      IF(I .EQ. 7) GO TO 1

      DO 12 J = 0,I3
   12 AC(J)=C(J)

    1 CONTINUE

C    REQUIRED ACCURACY NOT OBTAINED

      NC=64
      DELTA=0
      DO 13 J = 60,NC
   13 DELTA=DELTA+ABS(C(J))
      C(0)=HF*C(0)
      CALL MTLPRT(NAME,'E406.1','REQUIRED ACCURACY NOT OBTAINED')
      RETURN

C    REQUIRED ACCURACY OBTAINED
C    SUM NEGLECTED TERMS IN EXPANSION

    9 DELTA=0
      DO 7 K = J,I3
    7 DELTA=DELTA+CC(K)

C    CHECK IF FURTHER REDUCTION OF COEFFICIENTS IS POSSIBLE.

      NC=J-1
      REST=EPS-DELTA
      IF(REST .GT. 0) THEN
    8  IF(CC(NC) .GE. REST) GO TO 14
       DELTA=DELTA+CC(NC)
       REST=REST-CC(NC)
       NC=NC-1
       GO TO 8
      ENDIF
   14 C(0)=HF*C(0)
      RETURN
      END
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	*
	10	* e406cod.inc
	11	*
	12	DIMENSION C(0:128),CC(0:128),AC(0:64),BC(0:64)
	13
	14	PARAMETER (Z1 = 1, Z2 = 2, HF = Z1/2)
	15
	16	ALFA=HF*(B-A)
	17	BETA=HF*(B+A)
	18	C1=F(A)
	19	C2=F(B)
	20	AC(0)=C2+C1
	21	AC(1)=C2-C1
	22
	23	DO 1 I = 1,7
	24	I1=2**(I-1)
	25	I2=I1-1
	26	I3=2*I1
	27	C1=Z2/I1
	28	C2=PI/I1
	29
	30	DO 2 J = 0,I2
	31	2 C(J)=F(ALFACOS((J+HF)C2)+BETA)
	32
	33	C COMPUTE B-COEFFICIENTS
	34
	35	DO 3 J = 0,I2
	36	F1=J*C2
	37	F2=-HF*F1
	38
	39	C3=2*COS(F1)
	40	A2=0
	41	A1=0
	42	A0=C(I2)
	43	DO 4 K = I2-1,0,-1
	44	A2=A1
	45	A1=A0
	46	4 A0=C(K)+C3*A1-A2
	47	3 BC(J)=C1(A0COS(F1+F2)-A1*COS(F2))
	48	BC(I1)=0
	49
	50	C COMPUTE NEW C-COEFFICIENTS
	51
	52	DO 5 J = 0,I1
	53	C(J)=HF*(AC(J)+BC(J))
	54	5 C(I1+J)=HF*(AC(I1-J)-BC(I1-J))
	55
	56	C CHECK IF CALCULATION IS FINISHED
	57
	58	CMX=0
	59	DO 6 J = 0,I3
	60	CC(J)=ABS(C(J))
	61	6 CMX=MAX(CMX,CC(J))
	62
	63	IF(CMX .GT. 0) THEN
	64	CMX=1/CMX
65	CC(I3)=HF*CC(I3)
66	A0=CC(I2)*CMX
67	A1=CC(I1)*CMX
68	DO 10 J = I1+2,I3
69	A2=CC(J)*CMX
70	IF(A0 .LE. EPS .AND. A1. LE. EPS .AND. A2 .LE. EPS) GO TO 9
71
72	A0=A1
73	10 A1=A2
74	ENDIF
75
76	C DOUBLE THE NUMBER OF COEFFICIENTS.
77
78	IF(I .EQ. 7) GO TO 1
79
80	DO 12 J = 0,I3
81	12 AC(J)=C(J)
82
83	1 CONTINUE
84
85	C REQUIRED ACCURACY NOT OBTAINED
86
87	NC=64
88	DELTA=0
89	DO 13 J = 60,NC
90	13 DELTA=DELTA+ABS(C(J))
91	C(0)=HF*C(0)
92	CALL MTLPRT(NAME,'E406.1','REQUIRED ACCURACY NOT OBTAINED')
93	RETURN
94
95	C REQUIRED ACCURACY OBTAINED
96	C SUM NEGLECTED TERMS IN EXPANSION
97
98	9 DELTA=0
99	DO 7 K = J,I3
100	7 DELTA=DELTA+CC(K)
101
102	C CHECK IF FURTHER REDUCTION OF COEFFICIENTS IS POSSIBLE.
103
104	NC=J-1
105	REST=EPS-DELTA
106	IF(REST .GT. 0) THEN
107	8 IF(CC(NC) .GE. REST) GO TO 14
108	DELTA=DELTA+CC(NC)
109	REST=REST-CC(NC)
110	NC=NC-1
111	GO TO 8
112	ENDIF
113	14 C(0)=HF*C(0)
114	RETURN
115	END