C====================================================
C    PROGRAM FOR CALCULATION AMPLITUDE FACTOR
C    for endohedral atoms in one shell------
C    AUTHORS:A.S.Baltenkov,L.V.CHERNYSHEVA,M.YA.AMUSIA
C                VARIABLES
C      NES- NUMBER OF EXCITED STATES
C      NE5= NE+15, NE - NUMBER OF POINTS IN RADIAL MESH
C      S - NUMBER OF SHELLS
C      NW1 - MAX NUMBER OF FUNCTIONS IN FILE
C      RF-fullerene radius 
C      RF=6.720521 (C60), RF=13.5 (C240)
C      V0=-dl/2 binding potential
C      V0=0.4425 (C60) V0=0.5293 (C240)
C     dL=-2V0=-bet*(1+coth(bet*R)) (in atomic units). 
C     R fullerene radius, bet=sqrt(2*I), 
C     I=2.26 eV is affinity energy of the negative C60 ion.
C     ============================================
       SUBROUTINE FFACTOR1(NE,AS,LNO,NW1,M3,ALFA,BET,RF,V0,
     * AM,PH,EELRY,R1,R2,R8)
C =================================================


      INTEGER NE,NW1,KK,KK1,AS,LA,M3,LNO,NE5
      REAL*8 BET,ALFA,RF,PI,DL,V0,
     *       AA,BB,CC,AA1,BB1,CC1,DIS,DIS1,AJ,
     *       W1,W2,VN,VR,TF,KAP,
     *       EELRY(NW1),AM(NW1),PH(NW1), 
     *       R1(NE+15),R2(NE+15),R8(NE+15,*)


       PI=3.141592653
       NE5=NE+15
C-----------------------------------------
c      to calculate amplitude factor for endohedral atoms
c      in Dirac-bubble model (Jerusalem 2006)
C       RF=6.720521
C=====  RF is exp.fullerene radius in a.u.======
C       DL=-0.885
        DL=-2.*V0 
C====  DL is logarithmic derivative=====
C------------------------------------------
       DO 20 J=1,M3
       LA=AS+J
C       WRITE(6,52)LA,M3
  52   FORMAT(' LA,M3 ',2I5) 
      CALL DISKO1(LNO,J,2,NE5,R8(1,LA))
       CALL COLUMN(NE5,R8(1,LA),R2)
       LK=IDINT(R2(2))
C       WRITE(6,53)R2(1),LK
C  53   FORMAT(' E,L ',E15.8, I5) 
       IF(R2(1).GT.0.) THEN 
       KAP=DSQRT(R2(1))
       ELSE
       AM(J)=1.
       PH(J)=0.
       EELRY(J)=R2(1)
       GOTO  100
       END IF
       DO 51 II=3,NE+2
       AA=R1(II)/(ALFA*R1(II)+BET)
C       WRITE(6,52)AA
C  52   FORMAT(' AA ',E15.8) 
       R2(II)=R2(II)*DSQRT(AA)*DSQRT(PI*KAP)
  51   CONTINUE

C       WRITE(6,21) I,J,LK,KAP
C       WRITE(*,21) I,J,LK,KAP
  21   FORMAT('I,J, LK, KAPPA ',3I3,E15.8)
C       STOP
C       DO 501 IJ=1,NW1
C       WRITE(6,500)EPHRY(IJ)
C  500  FORMAT('EPHRy',E15.8)
C  501  CONTINUE
       DO 55 KK=3,NE+2
       IF(R1(KK)-RF)55,55,15
55     CONTINUE
C--- Parabolic-approximation ##########################
15     KK1=kK
       AA=(R2(KK1)-R2(KK1+1))/(R1(KK1+1)-R1(KK1-1))
     *                       /(R1(KK1)-R1(KK1+1))-
     *    (R2(KK1-1)-R2(KK1))/(R1(KK1+1)-R1(KK1-1))
     *                       /(R1(KK1-1)-R1(KK1))
       BB=(R2(KK1)-R2(KK1+1))/(R1(KK1)-R1(KK1+1))
     *                    -AA*(R1(KK1)+R1(KK1+1))
       CC=R2(KK1+1)-AA*R1(KK1+1)**2-BB*R1(KK1+1)
       KK1=kK+1
       AA1=(R2(KK1)-R2(KK1+1))/(R1(KK1+1)-R1(KK1-1))
     *                        /(R1(KK1)-R1(KK1+1))-
     *     (R2(KK1-1)-R2(KK1))/(R1(KK1+1)-R1(KK1-1))
     *                        /(R1(KK1-1)-R1(KK1))
       BB1=(R2(KK1)-R2(KK1+1))/(R1(KK1)-R1(KK1+1))
     *                    -AA1*(R1(KK1)+R1(KK1+1))
       CC1=R2(KK1+1)-AA1*R1(KK1+1)**2.-BB1*R1(KK1+1)
       AA=(AA+AA1)/2.
       BB=(BB+BB1)/2.
       CC=(CC+CC1)/2.
       DIS=BB*BB-4.*AA*CC
       DIS1=(2.*AA*Rf+BB)**2
       VR=AA*RF*RF+BB*RF+CC 
C       WRITE(6,11)DIS
11       FORMAT('Dis=',E15.8)      
C       IF(DIS)6,11,16
C 6      aJ=2./DSQRT(-dis)*DATAN((2.*AA*Rf+BB)/DSQRT(-dis))
C        GOTO 679
C 11     aJ=-2./(2.*AA*Rf+BB)
C        GOTO 679 
         IF(DIS.LT.0.)THEN
         aJ=2./DSQRT(-dis)*DATAN((2.*AA*Rf+BB)/DSQRT(-dis))
         END IF
         IF(DIS.EQ.0.)THEN
         aJ=-2./(2.*AA*Rf+BB)
         END IF
         IF(DIS.GT.0.)THEN
C        w1=2.*AA*Rf+BB-DSQRT(dis)
        w1=2.*AA*Rf+BB-DSQRT(dis)
        w2=2.*AA*Rf+BB+DSQRT(dis)
        aJ=1./DSQRT(dis)*DLOG(Dabs(w1/w2))
        END IF
       vn=2.*AA*Rf+BB+2.*AA*VR*AJ
       vn=-vn/dis
       vn=vn*KAP
C       WRITE(6,101) kap,rf,vr,vn
 101   FORMAT(' KAP rf,vr,vn ',4E15.8)
c--------Here vr and vn are regular and irregular 
C--------solutions at point r=Rf
       TF=VR**2/(VR*VN-KAP/DL)
       PH(J)=DACOS(1./DSQRT(1.+TF**2))
       AM(J)=(DCOS(PH(J))-DSIN(PH(J))*VN/VR)
       EELRY(J)=KAP*KAP
100     WRITE(6,95)J,EELRY(J),AM(J),PH(J)
       WRITE(*,95)J,EELRY(J),AM(J),PH(J)
 95     FORMAT(I3,3E15.8)
  20      CONTINUE
C---------------------------------------------
 600  RETURN
      END