Esope source files

   1 : C FARCA     SOURCE    PV        11/03/07    21:16:52     6885           
   2 :       FUNCTION FARCA(X,A,B,C,D,E,R,COST,SINT)
   3 :       IMPLICIT INTEGER(I-N)
   4 :       REAL*8 X,A,B,C,D,E
   5 :       REAL*8 T,ARG,DE
   6 :       REAL*8 RE,IM,SA,FARCA,AL,BE,DELTA
   7 :       REAL*8 C0,C1,C2,C3,C4,C5,C6,K1,K2
   8 :       REAL*8 R,COST,SINT,RES,MU,TD
   9 :       REAL*8 A2,B2,XI
  10 :       REAL*8 PREC
  11 : CCCCC CASTEM 2000 interdit COMPLEX*16.....
  12 : C      COMPLEX*16 Z1,FNUM1,DENRED,TMP
  13 :       COMPLEX*16 Z1,FNUM1,DENRED,TMP
  14 : C
  15 :       DATA PREC/1.D-8/
  16 : C
  17 :       SA=SQRT(A)
  18 :       AL=1.-D/SA
  19 :       BE=1.+D/SA
  20 :       DELTA=B**2-4.*A*C
  21 :       RE=-E+R*COST/(D-SA)
  22 :       IM=R*SINT/(D-SA)
  23 :       Z1=CMPLX(RE,IM)
  24 : C
  25 :       T=-SA*X+SQRT(A*X**2+B*X+C)
  26 :       DE=B-2.*SA*T
  27 :       C0=1./(4.*A)
  28 : C
  29 :       IF(ABS(IM).GT.PREC) THEN
  30 :        IF((ABS(DE).GT.PREC).AND.(ABS(DELTA).GT.PREC)) THEN
  31 :         ARG=(T-RE)/IM
  32 :         C1=B/(8.*A**1.5)-(B-2.*E*SA)/(4.*A*(D-SA))
  33 :         C4=(B**2-4.*A*C)**2/(8.*A**1.5)
  34 :         DENRED=(B-2.*SA*Z1)**3
  35 :         TMP=(Z1**2-C)**2
  36 :      &     *(2.*A*AL*Z1**2-2.*B*SA*AL*Z1+B**2-2.*A*C*BE)
  37 :         FNUM1=TMP/DENRED
  38 :         C5=IMAG(FNUM1)/IM
  39 :         C6=REAL(FNUM1)-C5*RE
  40 :         K1=-2.*C1*DELTA
  41 :         K1=K1-C4*(2.*B**3+24.*A*B*C)/DELTA**2
  42 :         K1=K1-2.*(C5*C*(B-2.*E*SA)-C6*(E*B-2.*C*SA))/R**2
  43 :         K2=-C0*DELTA/(2.*SA)
  44 :         K2=K2-C4*(3.*B**2+4.*A*C)/DELTA**2
  45 :         K2=K2+(C5*(E*B-2.*C*SA)-C6*(B-2.*E*SA))/(2.*SA*R**2)
  46 :         C3=B*K2-0.5*K1
  47 :         C2=K2-2.*B*C3/DELTA
  48 :         RES=0.5*C0*T**2+C1*T-C2*LOG(ABS(DE))/(2.*SA)
  49 :      &     +C3/(2.*SA*DE)+C4/(4.*SA*DE**2)
  50 :      &     +(C5*0.5*LOG(ABS(T**2-2.*RE*T+(RE**2+IM**2)))
  51 :      &     +(C6+RE*C5)/IM*ATAN(ARG))/(D-SA)
  52 :        ELSE
  53 :         A2=D+SA
  54 :         B2=B/(2.*SA)+E
  55 :         XI=A2*X+B2
  56 :         RES=(0.5*XI*(XI-4.*B2)+B2**2*LOG(ABS(XI)))/(A2**2)
  57 :        ENDIF
  58 :       ELSE
  59 :        IF((ABS(DE).GT.PREC).AND.(ABS(DELTA).GT.PREC)) THEN
  60 :         TD=(E*SA-0.5*B)/(D-SA)
  61 :         MU=-D*R-E
  62 :         C1=(3.*B-4.*MU/(D-SA))/(8.*A*SA)
  63 :         C5=2.*R**2
  64 :         C4=2.*MU**4/(A*SA)
  65 :         K1=-C1*DELTA-C4*(B**3+12.*A*B*C)/DELTA**2
  66 :         K1=K1+C5*TD*DELTA/(TD**2-C)
  67 :         K2=-C0*0.5*DELTA/SA-C4*(3.*B**2+4.*A*C)/DELTA**2
  68 :         K2=K2+C5*DELTA*0.5/((TD**2-C)*SA)
  69 :         C3=B*K2-K1
  70 :         C2=K2-2.*B*C3/DELTA
  71 :         RES=0.5*C0*T**2+C1*T-C2*LOG(ABS(DE))/(2.*SA)
  72 :      &     +C3/(2.*SA*DE)+C4/(4.*SA*DE**2)
  73 :      &     +C5*LOG(ABS(T-TD))
  74 :        ELSE
  75 :         A2=D+SA
  76 :         B2=B/(2.*SA)+E
  77 :         XI=A2*X+B2
  78 :         RES=(0.5*XI*(XI-4.*B2)+B2**2*LOG(ABS(XI)))/(A2**2)
  79 :        ENDIF
  80 :       ENDIF
  81 :       FARCA=0.5*RES
  82 :       RETURN
  83 :       END
  84 :  
  85 :  
  86 :