Esope source files

   1 : C FSEGA     SOURCE    CHAT      05/01/13    00:10:09     5004           
   2 :       FUNCTION FSEGA(GX,Y,X,A,D,E)
   3 :       IMPLICIT INTEGER(I-N)
   4 :       IMPLICIT REAL*8(A-H,O-Z)
   5 : *      REAL*8 GX,Y,X,A,D,E,aaa
   6 : *      REAL*8 SA,FSEGA
   7 :       REAL*8 MU1,MU2,MU3,MU4
   8 : *,T1,T2,T3,T4
   9 :       REAL*8 NORM2,NORM
  10 : *  ,ARG1,ALPHA1,SIG1,RAD
  11 : *      REAL*8 ARG2,ALPHA2,PRECIS
  12 : C
  13 :       PRECIS=1.D-12
  14 :       SA=SQRT(A)
  15 :       MU1=(D*GX+E)/SA
  16 :       MU2=D*X+E
  17 :       MU3=(GX-A*X-D*E)/SA
  18 :       MU4=GX-X
  19 : C
  20 :       SIG1=SIGN(1.D0,Y)
  21 :       NORM2=MU3**2+MU1**2
  22 :       NORM=SQRT(NORM2)
  23 :       RAD=SQRT(Y**2+NORM2)
  24 : C
  25 : 
  26 :       test2 = (abs(rad + mu2 ))/ ( abs(rad) + abs(mu2))
  27 :       test3= (abs( rad + mu3))/( abs(rad) + abs(mu3))
  28 :       aaa1 = rad + mu3
  29 :       aaa2 = rad + mu2
  30 :       if( test2.gt.1.e-13 .and.test3.gt.1.e-13.and.
  31 :      $  norm.ne.0.)then
  32 :         ARG1=(-MU3*RAD-NORM2)/(aaa1*NORM)
  33 :         ARG2=(-MU2*RAD-NORM2)/(aaa2*NORM)
  34 :          if( arg1.le.-1.d0) arg1=-1.d0
  35 :          if(arg1.gt.1.d0) arg1=1.d0
  36 :          if(arg2.le.-1.d0) arg2=-1.d0
  37 :          if(arg2.gt.1.d0) arg2=1.d0
  38 :         T1=MU1*(Y*LOG(RAD+MU3)-0.5D0*(Y-SIG1*ABS(MU1)*ASIN(ARG1)))
  39 :         T2=MU4*(Y*LOG(RAD+MU2)-0.5D0*(Y-SIG1*ABS(MU4)*ASIN(ARG2)))
  40 :         T3=0.5D0*(MU2*MU4+MU1*MU3)*SIG1*LOG((ABS(Y)+RAD)/NORM)
  41 :         IF(Y.NE.0.D0) THEN
  42 :           T4=-Y**2*ATAN((D-SA)*(SA*(MU4+MU1)+RAD)/Y)
  43 :         ELSE
  44 :           T4=0.D0
  45 :         ENDIF
  46 :         ALPHA1=-MU1*ABS(MU1)*ASIN(-1.)
  47 :         ALPHA2=-MU4*ABS(MU4)*ASIN(-1.)
  48 :         FSEGA=T1+T2+T3+T4+SIG1*0.5D0*(ALPHA2+ALPHA1)
  49 :       ELSE
  50 :         FSEGA=-Y**2*ATAN((D-SA)*SIG1)
  51 :       ENDIF
  52 :       RETURN
  53 :       END
  54 :  
  55 :  
  56 :  
  57 :