Sources Esope | Cast3M

fsega
C FSEGA     SOURCE    CHAT      05/01/13    00:10:09     5004
      FUNCTION FSEGA(GX,Y,X,A,D,E)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
*      REAL*8 GX,Y,X,A,D,E,aaa
*      REAL*8 SA,FSEGA
      REAL*8 MU1,MU2,MU3,MU4
*,T1,T2,T3,T4
      REAL*8 NORM2,NORM
*  ,ARG1,ALPHA1,SIG1,RAD
*      REAL*8 ARG2,ALPHA2,PRECIS
C
      PRECIS=1.D-12
      SA=SQRT(A)
      MU1=(D*GX+E)/SA
      MU2=D*X+E
      MU3=(GX-A*X-D*E)/SA
      MU4=GX-X
C
      SIG1=SIGN(1.D0,Y)
      NORM2=MU3**2+MU1**2
      NORM=SQRT(NORM2)
      RAD=SQRT(Y**2+NORM2)
C
 
      test2 = (abs(rad + mu2 ))/ ( abs(rad) + abs(mu2))
      test3= (abs( rad + mu3))/( abs(rad) + abs(mu3))
      aaa1 = rad + mu3
      aaa2 = rad + mu2
      if( test2.gt.1.e-13 .and.test3.gt.1.e-13.and.
     $  norm.ne.0.)then
        ARG1=(-MU3*RAD-NORM2)/(aaa1*NORM)
        ARG2=(-MU2*RAD-NORM2)/(aaa2*NORM)
         if( arg1.le.-1.d0) arg1=-1.d0
         if(arg1.gt.1.d0) arg1=1.d0
         if(arg2.le.-1.d0) arg2=-1.d0
         if(arg2.gt.1.d0) arg2=1.d0
        T1=MU1*(Y*LOG(RAD+MU3)-0.5D0*(Y-SIG1*ABS(MU1)*ASIN(ARG1)))
        T2=MU4*(Y*LOG(RAD+MU2)-0.5D0*(Y-SIG1*ABS(MU4)*ASIN(ARG2)))
        T3=0.5D0*(MU2*MU4+MU1*MU3)*SIG1*LOG((ABS(Y)+RAD)/NORM)
        IF(Y.NE.0.D0) THEN
          T4=-Y**2*ATAN((D-SA)*(SA*(MU4+MU1)+RAD)/Y)
        ELSE
          T4=0.D0
        ENDIF
        ALPHA1=-MU1*ABS(MU1)*ASIN(-1.)
        ALPHA2=-MU4*ABS(MU4)*ASIN(-1.)
        FSEGA=T1+T2+T3+T4+SIG1*0.5D0*(ALPHA2+ALPHA1)
      ELSE
        FSEGA=-Y**2*ATAN((D-SA)*SIG1)
      ENDIF
      RETURN
      END