C====MAIN PROGRAM ------------------
C---quadratic interpolation
C  for calcultion polarization factor G1(omega)
C    A@CN with one fullerene shell
C-------------------------------
      INTEGER   NE1,NE2
       REAL*8  RF,CO,EPH,KVINT,EPHRYP(100),
     *     EPHEVP(100),REGCN(100),IMGCN(100),
     *     REALFA(100),IMALFA(100),
     *     GCN(300),FCN(300),EPHRY1(300),
     *     GCN1(300),FCN1(300)


       NAMELIST/QPOL/ NE1,NE2,RF,EPHEVP,REALFA,IMALFA,
     *      EPHRY1

      OPEN(5,FILE='SYS$INPU',STATUS='OLD')
      READ(5,QPOL)
      OPEN(6,FILE='SYS$OUTP',STATUS='NEW')
C------- RF- theor.radius of fullerene---------
C      RF=6.720521
      CO=RF**3
C----GCN and FCN for plasm excitation-------
      WRITE(6, 1)RF
  1   FORMAT('              ENDOHEDRAL ATOMS '/
     *       '             POLARIZATION FACTOR  '/
     *       '     RADIUS OF FULLERENE  ',E15.8/
     *  4X,'OMEGA,eV',10X,'OMRGA,RY',10X,'ReG',
     *  10X,'IMG',10X,'G1_CN')
      DO 2 I=1,NE1
      EPHRYP(I)=EPHEVP(I)/13.6057
      REGCN(I)=1.-REALFA(I)/CO
      IMGCN(I)=-IMALFA(I)/CO
      GCN(I)=SQRT(REGCN(I)**2+IMGCN(I)**2)
      FCN(I)=ATAN2(IMGCN(I),REGCN(I))
      WRITE(6,3)EPHEVP(I),EPHRYP(I),REGCN(I),
     *          IMGCN(I),GCN(I)
  3   FORMAT(5E15.8)
  2   CONTINUE

        WRITE(6,5)
 5      FORMAT(' No    EPHRy',8X,'EPHeV',7X,'G1_CN',
     *             7X,'F1_CN')   
	DO 10 J=1,NE2
	EPH=EPHRY1(J)
        IF(EPH.GT.EPHRYP(1).AND.EPH.LT.EPHRYP(NE1))THEN
        GCN1(J)=KVINT(EPH,EPHRYP,GCN,1,NE1)
        FCN1(J)=KVINT(EPH,EPHRYP,FCN,1,NE1)
        EPH=EPH*13.5057
        WRITE(6,20)J,EPHRY1(J),EPH,GCN1(J),FCN1(J)
 20     FORMAT(I4,4E12.5)
        END IF
 10    CONTINUE
      STOP
      END
