C====MAIN PROGRAM ------------------
C---quadratic interpolation
C  for calcultion G12(omega) for 2 fullerene shells
C    A@C60C240
C-------------------------------
      INTEGER   NE1,NE2
       REAL*8  RF,C0,C1,C2,EPH,KVINT,EPHRYP(300),
     *   EPHEVP(300),REG12(300),IMG12(300),
     *   REALFA(300),IMALFA(300),
     *   G12(300),F12(300),EPHRY1(300),
     *   G121(300),F121(300)


       NAMELIST/QPOL2/ NE1,NE2,EPHEVP,REALFA,IMALFA,
     *      EPHRY1

      OPEN(5,FILE='SYS$INPU',STATUS='OLD')
      READ(5,QPOL2)
      OPEN(6,FILE='SYS$OUTP',STATUS='NEW')
      WRITE(6, 1)

  1   FORMAT('              ENDOHEDRAL ATOMS '/
     *     '  POLARIZATION FACTOR FOR 2 FULLERENE SHELLS')
      WRITE(6, 105)
105   FORMAT(4X,'OMEGA,eV',10X,'Realfa',10X,'Imalfa')
      DO 101 I=1,NE1 
      WRITE(6,100)EPHEVP(I),REALFA(I),IMALFA(I)
100   FORMAT(3E15.8)
101   CONTINUE
C      C0=2**4.
C      WRITE(6,101)C0
C101   FORMAT(' C0**5',E15.8)
C------- 6.62- theor.radius of fullerene---------
      RF=6.720521
      C0=1.51
      C1=0.38
      C2=0.0635
C----G12 and F12 for plasm excitation-------
      WRITE(6, 102)
102   FORMAT(4X,'OMEGA,eV',3X,'OMEGA,RY',6X,'ReGC60C240',
     *  6X,'IMGC60C240',6X,'G_C60C240',6X,'phase')
      DO 2 I=1,NE1
      EPHRYP(I)=EPHEVP(I)/13.6057
      REG12(I)=(RF**12.-C0*RF**9.*REALFA(I)+
     *     (C0*C1-2.*C2)*RF**6.*(REALFA(I)**2-IMALFA(I)**2)+
     *   C0*C2*RF**3.*REALFA(I)*(REALFA(I)**2+IMALFA(I)**2)+
     *   C2*(C2-C0*C1)*(REALFA(I)**2+IMALFA(I)**2)**2)/
     *  (RF**12.-2.*c2*RF**6.*(REALFA(I)**2-IMALFA(I)**2)+
     *  C2**2*(REALFA(I)**2+IMALFA(I)**2)**2)
      IMG12(I)=C0*RF**3.*IMALFA(I)/
     *   (RF**12.-2.*C2*RF**6.*(REALFA(I)**2-IMALFA(I)**2)+
     *   C2**2*(REALFA(I)**2+IMALFA(I)**2)**2)*
     *   (RF**6.-2.*C1*RF**3.*REALFA(I)+
     *   C2*(REALFA(I)**2+IMALFA(I)**2))
      G12(I)=SQRT(REG12(I)**2+IMG12(I)**2)
      F12(I)=ATAN2(IMG12(I),REG12(I))
      WRITE(6,3)EPHEVP(I),EPHRYP(I),REG12(I),
     *          IMG12(I),G12(I),F12(I)
  3   FORMAT(E11.4,5E15.8)
  2   CONTINUE

        WRITE(6,5)
 5      FORMAT(' No    EPHRy        EPHeV ',
     *   '   GC60C240   FC60C240  ' )    
	DO 10 J=1,NE2
	EPH=EPHRY1(J)
        IF(EPH.GT.EPHRYP(1).AND.EPH.LT.EPHRYP(NE1))THEN
        G121(J)=KVINT(EPH,EPHRYP,G12,1,NE1)
        F121(J)=KVINT(EPH,EPHRYP,F12,1,NE1)
        EPH=EPH*13.6057
        WRITE(6,20)J,EPHRY1(J),EPH,G121(J),F121(J)
 20     FORMAT(I4,4E12.5)
        END IF
 10    CONTINUE
      STOP
      END
