PROGRAM NEW6_ADAPT
      IMPLICIT none
      INTEGER*4 N,M,NAT,nw,inw,iq,nadapt,iwr
      real*8 EXCA,wmin,wmax,EXCNM,EXCCM,EXCAU,CM,wd,radapt
      real*8,allocatable::S(:,:),ALPHA(:,:,:),A(:,:,:,:),
     1G(:,:,:),R(:,:),alphav(:,:,:),E(:),
     1Gi(:,:,:),ALPHAi(:,:,:),Ai(:,:,:,:),we(:)
      LOGICAL LPOLAR,lpro,lgpro,lapro,
     1luseA,luseG,lor,lintA,lintG,lcompl,ladapt
C
      CM=219470.0d0
      WRITE(6,6666)
6666  FORMAT(/,/,
     A'          * * * * * * * * * * * * * * * * * * * * * *',/,
     1'        *                                         *  ',/,
     2'      *           RAMAN OPTICAL ACTIVITY        *    ',/,
     3'        *                                         *  ',/,
     4'          * * * * * * * * * * * * * * * * * * * * * *',/,
     8'                                                     ',/,
     8'            Normal mode - adaptive exc. frequencies  ',/,
     8'                                                     ',/,
     8'                       Petr Bour 1997-2025           ',/,/,/,
     A'INPUT FILES: ROA.OPT or NEW6.OPT,FILE.TTT,F.INP',/,
     A'   For adaptive frequencies:',/,
     A'   TTT.LST   list of the files with polarizability dervs.',/,/,
     B'OUTPUT     : ROA.TAB ',/)

      CALL IOPAR(LPOLAR,EXCA,wmin,wmax,lpro,lgpro,lapro,
     1luseA,luseG,lor,lintA,lintG,lcompl,ladapt,radapt,nadapt,iwr)
      EXCNM=EXCA/10.0d0
      EXCCM=1.0d7/EXCNM
      EXCAU=EXCCM/CM

      open(9,file='F.INP')
      read(9,*)M,N,NAT
      close(9)
      allocate(S(N,N),R(3,NAT),E(N),ALPHA(N,3,3),A(N,3,3,3),G(N,3,3),
     1alphav(N,3,3),ALPHAi(N,3,3),Ai(N,3,3,3),Gi(N,3,3))
      CALL READS(N,S,R,E,M)

      if(ladapt)then
c      determine number of frequencies
       nw=inw('TTT.LST')
       allocate(we(nw))
c      fill available excitation frequencies:
       call fille(we,nw)
       call openfiles(lor)
c      for each normal mode
       do 1 iq=1,M
c      desired exc. frequency for this mode:
       wd=EXCAU-radapt*E(iq)
c      make the tensors for it:
       call mkt(N,nw,wd,we,ALPHA,A,G,ALPHAi,Ai,Gi,lcompl,iq,E(iq),
     1 nadapt,iwr)
       CALL TRTEN(S,N,M,ALPHA,A,G)
       CALL TRTEN(S,N,M,ALPHAi,Ai,Gi)
       if(lcompl)then
        CALL RSi(N,ALPHA,ALPHAi,G,Gi,A,Ai,EXCA,wmin,wmax,E,
     1  luseA,luseG,lor,iq)
       else
        CALL RS(N,ALPHA,G,A,EXCA,wmin,wmax,E,
     1  luseA,luseG,lor,iq)
       endif
1      continue
       call closefiles(lor)
      else
       CALL READTEN(N,NAT,ALPHA,G,A,0)
       if(lcompl)then
        CALL READTEN(N,NAT,ALPHAi,Gi,Ai,3)
        WRITE(6,*)'(Complex part)'
       endif
C
       IF(LPOLAR)THEN
c       set A, G to zero:
        CALL TTTT(N,ALPHA,A,G,NAT,3,R)
c       make common A, G from alpha:
        CALL TTTT(N,ALPHA,A,G,NAT,2,R)
        write(6,*)' A, G made from alpha !!!'
       ENDIF

       if(.not.lintA.or..not.lintG)then
c       set internal A, G to zero:
c       1) transform A, G to atomic origins:
        CALL TTTT(N,ALPHA,A,G,NAT,1,R)
c       2) delete them
        if(.not.lintG)CALL TTTT(N,ALPHA,A,G,NAT,4,R)
        if(.not.lintA)CALL TTTT(N,ALPHA,A,G,NAT,5,R)
c       3) transform A, G back to common origin:
        CALL TTTT(N,ALPHA,A,G,NAT,2,R)
       endif
c
       call rules(ALPHA,A,G,NAT)
C
C      transform the tensors into internal coordinate derivatives
       CALL TRTEN(S,N,M,ALPHA,A,G)
       if(lcompl)then
        CALL TRTEN(S,N,M,ALPHAi,Ai,Gi)
        alphav=ALPHAi
        CALL WTQi(N,M,ALPHA,A,G,ALPHAi,Ai,Gi,alphav,E,WMIN,WMAX)
       else
C       write used normal mode tensor derivatives into FILE.Q.TTT
        alphav=ALPHA
        CALL WRITETTTQ(N,M,ALPHA,A,G,alphav,E,WMIN,WMAX)
       endif
c
       if(lcompl)then
        CALL DORAMANi(N,M,ALPHA,ALPHAi,G,Gi,A,Ai,EXCA,wmin,wmax,E,
     1  luseA,luseG,lor)
       else
        CALL DORAMAN(N,M,ALPHA,G,A,EXCA,wmin,wmax,E,
     1  luseA,luseG,lor)
       endif
      endif
      
      END
C     ==========================================
      SUBROUTINE WTQi(N,M,ALPHA,A,G,ALPHAi,Ai,Gi,av,E,WMIN,WMAX)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION ALPHA(N,3,3),G(N,3,3),A(N,3,3,3),
     1ALPHAi(N,3,3),Gi(N,3,3),Ai(N,3,3,3),av(N,3,3),E(*)
c
      CM=219470.0d0
      im=0
      DO 4 II=1,M
4     if(E(II)*CM.GE.WMIN.AND.E(II).LE.WMAX)im=im+1
      open(22,file='FILE.Q.TTT')
      WRITE(22,2003)im
2003  FORMAT(' ROA tensors, normal modes derivatives',/,I4,' modes',/,
     1' The electric-dipolar electric-dipolar polarizability:',/,
     2        ' mode e(cm-1)   jx           jy           jz')
      im=0
      DO 1 II=1,M
      if(E(II)*CM.GE.WMIN.AND.E(II).LE.WMAX)then
       im=im+1
       write(22,221)im,E(II)*CM
221    format(i5,f12.2)
       DO 11 I=1,3
11     WRITE(22,2001)I,(ALPHA(II,I,J),J=1,3),(ALPHAi(II,I,J),J=1,3)
2001   FORMAT(I3,6E14.6)
      endif
1     continue

      WRITE(22,2004)
2004  FORMAT(' The electric dipole magnetic dipole polarizability:',/,
     1      ' mode  e(cm-1)  jx(Bx)       jy(By)       jz(Bz)')
      im=0
      DO 2 II=1,M
      if(E(II)*CM.GE.WMIN.AND.E(II).LE.WMAX)then
       im=im+1
       write(22,221)im,E(II)*CM
       DO 21 I=1,3
21     WRITE(22,2001)I,(G(II,I,J),J=1,3),(Gi(II,I,J),J=1,3)

c      last inner index magnetic
      endif
2     continue

      WRITE(22,2005)
2005  FORMAT(' The electric dipole electric quadrupole polarizability:',
     2/,     ' mode  e(cm-1)  kx           ky           kz')
      im=0
      DO 3 II=1,M
      if(E(II)*CM.GE.WMIN.AND.E(II).LE.WMAX)then
       im=im+1
       write(22,221)im,E(II)*CM
       DO 31 I=1,3
       DO 31 J=1,3
31     WRITE(22,2002)I,J,(A(II,I,J,K),K=1,3),(Ai(II,I,J,K),K=1,3)

2002   FORMAT(2I3,6F14.6)
      endif
3     continue

      WRITE(22,2009)
2009  FORMAT('The velocity form of alpha',/,
     1' mode e(cm-1)    vx           vy           vz')
      im=0
      DO 5 II=1,M
      if(E(II)*CM.GE.WMIN.AND.E(II).LE.WMAX)then
       im=im+1
       write(22,221)im,E(II)*CM
       DO 51 I=1,3
51     WRITE(22,2001)I,(av(II,I,J),J=1,3),(av(II,I,J),J=1,3)
      endif
5     continue

      close(22)
      RETURN
      END
C     ==========================================
      SUBROUTINE WRITETTTQ(N,M,ALPHA,A,G,alphav,E,WMIN,WMAX)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION ALPHA(N,3,3),G(N,3,3),A(N,3,3,3),
     1alphav(N,3,3),E(*)
c
      CM=219470.0d0
      im=0
      DO 4 II=1,M
4     if(E(II)*CM.GE.WMIN.AND.E(II).LE.WMAX)im=im+1
      open(22,file='FILE.Q.TTT')
      WRITE(22,2003)im
2003  FORMAT(' ROA tensors, normal modes derivatives',/,I4,' modes',/,
     1' The electric-dipolar electric-dipolar polarizability:',/,
     2        ' mode e(cm-1)   jx           jy           jz')
      im=0
      DO 1 II=1,M
      if(E(II)*CM.GE.WMIN.AND.E(II).LE.WMAX)then
       im=im+1
       write(22,221)im,E(II)*CM
221    format(i5,f12.2)
       DO 11 I=1,3
11     WRITE(22,2001)I,(ALPHA(II,I,J),J=1,3)
2001   FORMAT(I3,3F13.7)
      endif
1     continue

      WRITE(22,2004)
2004  FORMAT(' The electric dipole magnetic dipole polarizability:',/,
     1      ' mode  e(cm-1)  jx(Bx)       jy(By)       jz(Bz)')
      im=0
      DO 2 II=1,M
      if(E(II)*CM.GE.WMIN.AND.E(II).LE.WMAX)then
       im=im+1
       write(22,221)im,E(II)*CM
       DO 21 I=1,3
21     WRITE(22,2001)I,(G(II,I,J),J=1,3)
c      last inner index magnetic
      endif
2     continue

      WRITE(22,2005)
2005  FORMAT(' The electric dipole electric quadrupole polarizability:',
     2/,     ' mode  e(cm-1)  kx           ky           kz')
      im=0
      DO 3 II=1,M
      if(E(II)*CM.GE.WMIN.AND.E(II).LE.WMAX)then
       im=im+1
       write(22,221)im,E(II)*CM
       DO 31 I=1,3
       DO 31 J=1,3
31     WRITE(22,2002)I,J,(A(II,I,J,K),K=1,3)
2002   FORMAT(2I3,3F13.7)
      endif
3     continue

      WRITE(22,2009)
2009  FORMAT('The velocity form of alpha',/,
     1' mode e(cm-1)    vx           vy           vz')
      im=0
      DO 5 II=1,M
      if(E(II)*CM.GE.WMIN.AND.E(II).LE.WMAX)then
       im=im+1
       write(22,221)im,E(II)*CM
       DO 51 I=1,3
51     WRITE(22,2001)I,(ALPHAV(II,I,J),J=1,3)
      endif
5     continue

      close(22)
      RETURN
      END
C     ==========================================
      FUNCTION EPS(I,J,K)
      IMPLICIT INTEGER*4 (I-N)
      REAL*8 EPS
      EPS=0.0D0
      IF(I.EQ.1)then
       if(J.EQ.2.AND.K.EQ.3)EPS= 1.0D0
       IF(J.EQ.3.AND.K.EQ.2)EPS=-1.0D0
      endif
      IF(I.EQ.2)then
       IF(J.EQ.3.AND.K.EQ.1)EPS= 1.0D0
       IF(J.EQ.1.AND.K.EQ.3)EPS=-1.0D0
      endif
      IF(I.EQ.3)then
       IF(J.EQ.1.AND.K.EQ.2)EPS= 1.0D0
       IF(J.EQ.2.AND.K.EQ.1)EPS=-1.0D0
      endif
      RETURN
      END
C     ==============================================================
      SUBROUTINE READS(N,S,R,E,M)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION S(N,N),E(N),R(3,N/3)
      S=0.0d0
      CM=219470.0d0
      SFAC=0.0234280d0
C     1 a.u.= 2.1947E5 cm^-1
      OPEN(4,FILE='F.INP',STATUS='OLD')
      read(4,*)M,N,nat
      do 1 i=1,NAT
1     read(4,*)R(1,i),(R(ix,i),ix=1,3)
      read(4,*)
      DO 2 I=1,NAT
      DO 2 J=1,M
2     read(4,*)(s(3*(i-1)+ix,j),ix=1,2),(s(3*(i-1)+ix,j),ix=1,3)
      read(4,*)
      READ(4,4000)(E(I),I=1,M)
4000  FORMAT(6F11.6)
      close(4)
      E=E/CM
      S=S*SFAC
      RETURN
      end

      SUBROUTINE READTTT(N,NAT,P,G,A,ic,fn)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION P(N,3,3),G(N,3,3),A(N,3,3,3)
      character*(*) fn
      OPEN(2,FILE=fn,STATUS='OLD')
      READ(2,*)
      READ(2,*)NAT
      if(NAT.gt.N/3)call report('too many atoms')
      READ(2,*)
      READ(2,*)
      DO 1 I=1,3
      READ(2,*)
      DO 1 L=1,NAT
      DO 1 IX=1,3
      M=3*(L-1)+IX
1     READ(2,*)(P(M,I,J),J=1,2),(P(M,I,J),J=1,ic),(P(M,I,J),J=1,3)
      READ(2,*)
      READ(2,*)
      DO 2 I=1,3
      READ(2,*)
      DO 2 L=1,NAT
      DO 2 IX=1,3
      M=3*(L-1)+IX
c     last index magnetic, inner
2     READ(2,*)(G(M,I,J),J=1,2),(G(M,I,J),J=1,ic),(G(M,I,J),J=1,3)
      READ(2,*)
      READ(2,*)
      DO 3 I=1,3
      DO 3 J=1,3
      READ(2,*)
      DO 3 L=1,NAT
      DO 3 IX=1,3
      M=3*(L-1)+IX
3     READ(2,*)(A(M,I,J,K),K=1,2),(A(M,I,J,K),K=1,ic),
     1(A(M,I,J,K),K=1,3)
      CLOSE(2)
      RETURN
      END
C     ==============================================================
C
      SUBROUTINE READTEN(N,NAT,P,G,A,ic)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION P(N,3,3),G(N,3,3),A(N,3,3,3)
      OPEN(2,FILE='FILE.TTT',STATUS='OLD')
      READ(2,*)
      READ(2,*)NAT
      if(NAT.gt.N/3)call report('too many atoms')
      READ(2,*)
      READ(2,*)
      DO 1 I=1,3
      READ(2,*)
      DO 1 L=1,NAT
      DO 1 IX=1,3
      M=3*(L-1)+IX
1     READ(2,*)(P(M,I,J),J=1,2),(P(M,I,J),J=1,ic),(P(M,I,J),J=1,3)
      READ(2,*)
      READ(2,*)
      DO 2 I=1,3
      READ(2,*)
      DO 2 L=1,NAT
      DO 2 IX=1,3
      M=3*(L-1)+IX
c     last index magnetic, inner
2     READ(2,*)(G(M,I,J),J=1,2),(G(M,I,J),J=1,ic),(G(M,I,J),J=1,3)
      READ(2,*)
      READ(2,*)
      DO 3 I=1,3
      DO 3 J=1,3
      READ(2,*)
      DO 3 L=1,NAT
      DO 3 IX=1,3
      M=3*(L-1)+IX
3     READ(2,*)(A(M,I,J,K),K=1,2),(A(M,I,J,K),K=1,ic),
     1(A(M,I,J,K),K=1,3)
      CLOSE(2)
      WRITE(6,*)NAT,' atoms, tensors read in from FILE.TTT'
      RETURN
      END
C     ==============================================================
      SUBROUTINE IOPAR(LPOLAR,EXCA,wmin,wmax,lpro,lgpro,lapro,
     1luseA,luseG,lor,lintA,lintG,lcompl,ladapt,radapt,nadapt,iwr)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8(A-H,O-Z)
      LOGICAL LPOLAR,lpro,lex,lgpro,lapro,
     1luseA,luseG,lor,lintA,lintG,lcompl,ladapt
      CHARACTER*15 STRING,st
c
      wmin=0.0d0
      wmax=1.0d90
      lpro=.false.
      lapro=.false.
      lcompl=.false.
      lgpro=.false.
      LPOLAR=.FALSE.
c     amount if output:
      IWR=0
c     adaptive excitation frequency w -> w- radapt x wI :
      ladapt=.true.
c     the adaptation factor
      radapt=0.5d0
c     interpolate polynomila degree (0 closest neighbors):
      nadapt=0
c     excitation light wavelength in nm:
      EXCNM=532.0d0
c     kind of ROA:
c     0 - 180 deg SCP (ordinary)
c     1 - 90 deg polarized
c     2 - 90 deg depolarized
c     use A-tensor derivatives:
      luseA=.true.
c     use G'-tensor derivatives:
      luseG=.true.
c     use intrinsic part of A-tensor derivatives:
      lintA=.true.
c     use intrinsic part of G'-tensor derivatives:
      lintG=.true.
c     generate also spectra for oriented sampels:
      lor=.false.

      INQUIRE(FILE='ROA.OPT',exist=lex)
      if(lex)then
       WRITE(6,*)' FILE ROA.OPT FOUND ...'
       OPEN(2,FILE='ROA.OPT')
      else
       INQUIRE(FILE='NEW6.OPT',exist=lex)
       if(lex)then
        OPEN(2,FILE='NEW6.OPT')
        WRITE(6,*)' FILE NEW6.OPT FOUND ...'
       else
        goto 99
       endif
      endif
   
2001  READ(2,1111)STRING
1111  FORMAT(A15)
      IF(STRING(1:4).EQ.'ADAP')READ(2,*)ladapt
      IF(STRING(1:5).EQ.'POLAR')READ(2,*)LPOLAR
      IF(STRING(1:4).EQ.'WMIN')READ(2,*)wmin
      IF(STRING(1:4).EQ.'WMAX')READ(2,*)wmax
      IF(STRING(1:5).EQ.'EXCNM')READ(2,*)EXCNM
      IF(STRING(1:5).EQ.'COMPL')READ(2,*)lcompl
      IF(STRING(1:3).EQ.'PRO')READ(2,*)lpro
      IF(STRING(1:3).EQ.'IWR')READ(2,*)iwr
      IF(STRING(1:4).EQ.'GPRO')READ(2,*)lgpro
      IF(STRING(1:4).EQ.'APRO')READ(2,*)lapro
      IF(STRING(1:4).EQ.'INTA')READ(2,*)lintA
      IF(STRING(1:4).EQ.'INTG')READ(2,*)lintG
      IF(STRING(1:4).EQ.'USEA')READ(2,*)luseA
      IF(STRING(1:4).EQ.'USEG')READ(2,*)luseG
      IF(STRING(1:5).EQ.'ORIEN')READ(2,*)lor
      IF(STRING(1:5).EQ.'RADAP')READ(2,*)radapt
      IF(STRING(1:5).EQ.'NADAP')READ(2,*)nadapt
      IF(STRING(1:3).EQ.'END')GOTO 2002
      backspace 2
      read(2,1111)st
      if(st.eq.STRING)then
       write(6,1111)st
       call report('Unknown option')
      endif
      GOTO 2001
2002  CLOSE(2)
     
99    write(6,6009)EXCNM
 6009 format(' Excitation ',f5.1,' nm')
      EXCA=EXCNM*10.0D0
      RETURN
      END

      SUBROUTINE RS(N,ALPHA,G,A,EXCA,wmin,wmax,E,
     1luseA,luseG,lor,iq)
c
c     supposedly G/W is supplied
c
      IMPLICIT none
      integer*4 IQ,I,J,K,L,N
      real*8 ALPHA(N,3,3),G(N,3,3),A(N,3,3,3),E(*),EXCA,wmin,wmax,
     1CM,ECM,SAL0,SAL1,SAG0,SAG1,SA1,SPAL,SDCP180,SDCPII180,
     1S180,D180,D,S90X,D90X,DX,S90Z,D90Z,DZ,YDY,YDX,P1,down,doc,
     1ro,AMU,BOHR,SpA3,beta2,gpisvejc,alphag,betaa,
     1betag,alpha2,roa3,roa2,ram2,CID2,roa1,
     1ram1,EPS,CID1,IRR180,IRL180
      logical luseA,luseG,lor
c
      call init(N,CM,luseA,luseG,A,G)
c
      ECM=E(IQ)*CM
      IF(ABS(ECM).GT.0.01d0)THEN
      if(wmin.eq.0.0d0.or.ECM.gt.wmin)then
      if(wmax.eq.0.0d0.or.ECM.lt.wmax)then
       SAL0=0.0d0
       SAL1=0.0d0
       SAG0=0.0d0
       SAG1=0.0d0
       SA1=0.0d0
       SPAL=0.0d0
       DO 3 I=1,3
       SPAL=SPAL+ALPHA(IQ,I,I)
       DO 3 J=1,3
       SAL0=SAL0+ALPHA(IQ,I,I)*ALPHA(IQ,J,J)
       SAL1=SAL1+ALPHA(IQ,I,J)*ALPHA(IQ,I,J)
       SAG0=SAG0+ALPHA(IQ,I,I)*    G(IQ,J,J)
       SAG1=SAG1+ALPHA(IQ,I,J)*    G(IQ,I,J)
       DO 3 K=1,3
       DO 3 L=1,3
3      SA1=SA1+EPS(I,K,L)*ALPHA(IQ,I,J)*A(IQ,K,L,J)/3.0d0
c
c      IL+IR, backscattering, DCP_I,DCP_II
       SDCP180  =6.0d0*(6.0d0*SAL1-2.0d0*SAL0)
       SDCPII180=6.0d0*(      SAL1+3.0d0*SAL0)
c
c      IL+/-IR, backscattering
       S180=7.0d0*SAL1+SAL0
c      write(6,*)ECM,SAL0,SAL1,SAG0,SAG1,SA1
       D180=8.0d0*(3.0d0*SAG1-SAG0+SA1)
       D=D180
       IF(S180.gt.1.0d-9)D=D/S180
c
c      IL+/-IR, 90o, x
       S90X=7.0d0*SAL1+SAL0
       D90X=2.0d0*(7.0d0*SAG1+SAG0+SA1)
       DX=D90X
       IF(S90X.gt.1.0d-9)DX=DX/S90X
c
c      IL+/-IR, 90o, z
       S90Z=2.0d0*(3.0d0*SAL1-SAL0)
       D90Z=4.0d0*(3.0d0*SAG1-SAG0-SA1)
       DZ=D90Z
       IF(S90Z.gt.1.0d-9)DZ=DZ/S90Z
c
c      polarized backscattering components (YDY+YDX=S180):
       YDY=3.0d0*SAL1-SAL0
       YDX=4.0d0*SAL1+2.0d0*SAL0
       if(YDX.gt.1.0d-9)then
        P1=YDY/YDX
       else
        P1=0.0d0
       endif
      
c      calculate degree of circularity
       down=        7.0d0*SAL1+ 1.0d0*SAL0
       doc=5.0d0*(  1.0d0*SAL1- 1.0d0*SAL0)
       if(down.ne.0.0d0)doc=doc/down

c      calculate depolarization
       down= 6.0d0*SAL0-12.0d0*SAL1
       ro  =-3.0d0*SAL0+ 9.0d0*SAL1
       if(down.ne.0.0d0)ro=ro/down
c
       AMU=1822.0d0
       BOHR=0.529177d0
       Alpha2=SAL0/9.0d0*(AMU*BOHR**4)
       SpA3=SPAL/3.0d0*dsqrt((AMU*BOHR**4))
       beta2=0.5d0*(3*SAL1-SAL0)*(AMU*BOHR**4)
       SDCP180  =SDCP180  *(AMU*BOHR**4)
       SDCPII180=SDCPII180*(AMU*BOHR**4)
       YDY=6.0d0*YDY*(AMU*BOHR**4)
       YDX=6.0d0*YDX*(AMU*BOHR**4)
c      gpisvejc=980.0
       gpisvejc=(AMU*BOHR**5)*1.0d4*2.0d0*4.0d0*atan(1.0d0)/EXCA
c
c      Barron's tensor invariants:
c      alpha G = [(1/3)Sp alpha] [(1/3)Sp G'] etc
c      beta(A)^2=Delta2 in Gaussian:
       alphag=SAG0/9.0d0                *gpisvejc
       betaa=SA1*3.0d0/2.0d0            *gpisvejc
       betag=(3.0d0*SAG1-SAG0)/2.0d0    *gpisvejc
       WRITE(44,4305)IQ,ECM,Alpha2,beta2,alphag,betag,betaa
4305   FORMAT(i6,f10.2,5f11.3)
c
c      DCPI 180
c      ram3=24.0d0*beta2 = 12 (3al_ab al_ab -al_aa al_bb)
       roa3=8.0d0*(12.0d0*betag+4.0d0*betaa)
c
c      polarized ROA,90deg:
c      roa90x=45.0d0*alphag+7.0d0*betag+betaa
c      ram90x=0.5d0*(45.0d0*alpha2+7.0d0*beta2)
c      depolarized ROA,90deg:
c      roa90z=6.0d0*(betag-betaa/3.0d0)
c      ram90z=0.5d0*6.0d0*beta2
c
c      unknown:
       S90X=S90X*(AMU*BOHR**4)
       D90X=D90x*gpisvejc
       DX=D90X
       IF(S90X.gt.1.0d-9)DX=DX/S90X
c
c      ICP(90)=90Z
       roa2=8.0d0*(3.0d0*betag-betaa)
       ram2=12.0d0*beta2
       CID2=0.0d0
       if(ram2.gt.1.0d-9)CID2=roa2/ram2
c
c      ICP/SCP(180):
c      8/c*(12*beta(G')^2 + 4*beta(A)^2):
       roa1=96.0d0*(betag+betaa/3.0d0)
c      ram1 = 6(7 al_ab al_ab + al_aa al_bb)
       ram1=4.0d0*(45.0d0*SpA3**2+7.0d0*beta2)
       CID1=0.0d0
       if(ram1.gt.1.0d-9)CID1=roa1/ram1

       IRR180 = 24.0d0*beta2
       IRL180 =180.0d0*SpA3**2+4.0d0*beta2


       WRITE(4,3000)IQ,ECM, YDX,YDY,ram1, CID2,DX,CID1, P1,
     1 roa1,doc,SDCP180,roa3,SDCPII180,IRR180,IRL180
3000   FORMAT(I5,f9.2,6e12.4,f6.3,e12.4,f7.3,7e12.4)
       if(lor)then
        do 6 i=1,3
        call setjk(i,j,k)
c       S0(x)=(1/2)(azz azz + ayy ayy + azy azy + azy azy)
c       S0(y)=(1/2)(axx axx + azz azz + axz axz + axz axz)
c       S0(z)=(1/2)(axx axx + ayy ayy + axy axy + axy axy)
c                    jj  jj    kk  kk    jk  jk    jk  jk
        ram1=180.0d0*0.5d0*AMU*BOHR**4*(
     1  ALPHA(IQ,j,j)*ALPHA(IQ,j,j)+ALPHA(IQ,k,k)*ALPHA(IQ,k,k)+
     1  ALPHA(IQ,j,k)*ALPHA(IQ,j,k)+ALPHA(IQ,j,k)*ALPHA(IQ,j,k))
c      -2ayx(Gxy+Gyx)+(ayy-axx)(Gxx-Gyy)
c         kj  jk  kj    kk  jj   jj  kk
c       +(2w/3)(ayyAx zy-axxAy zx+ayxAx zx-axyAy zy)
c                kk j ik  jj k ij  kj j ij  jk k ik
        roa1=-180.0d0*gpisvejc*(
     1  -2.0d0*ALPHA(IQ,j,k)*(G(IQ,j,k)+G(IQ,k,j))
     1  +(ALPHA(IQ,k,k)-ALPHA(IQ,j,j))*(G(IQ,j,j)-G(IQ,k,k))
     1  +2.0d0/3.0d0*
     1  (ALPHA(IQ,k,k)*A(IQ,j,i,k)-ALPHA(IQ,j,j)*A(IQ,k,i,j)
     1  +ALPHA(IQ,k,j)*A(IQ,j,i,j)-ALPHA(IQ,j,k)*A(IQ,k,i,k)))
6       WRITE(9+i,3000)IQ,ECM,YDX,YDY,ram1,CID2,DX,CID1,P1,roa1
       endif
c
      endif
      endif
      ENDIF

      RETURN
      END


      SUBROUTINE DORAMAN(M,N,ALPHA,G,A,EXCA,wmin,wmax,E,
     1luseA,luseG,lor)
c
c     supposedly G/W is supplied
c
      IMPLICIT none
      integer*4 IQ,I,J,K,L,N,M
      real*8 ALPHA(N,3,3),G(N,3,3),A(N,3,3,3),E(*),EXCA,wmin,wmax,
     1CM,ECM,SAL0,SAL1,SAG0,SAG1,SA1,SPAL,SDCP180,SDCPII180,
     1S180,D180,D,S90X,D90X,DX,S90Z,D90Z,DZ,YDY,YDX,P1,down,doc,
     1ro,AMU,BOHR,SpA3,beta2,gpisvejc,alphag,betaa,
     1betag,alpha2,roa3,roa2,ram2,CID2,roa1,
     1ram1,EPS,CID1,IRR180,IRL180
      logical luseA,luseG,lor
c
      call init(N,CM,luseA,luseG,A,G)
      call openfiles(lor)
c
      DO 1 IQ=1,M
      ECM=E(IQ)*CM
      IF(ABS(ECM).GT.0.01d0)THEN
      if(wmin.eq.0.0d0.or.ECM.gt.wmin)then
      if(wmax.eq.0.0d0.or.ECM.lt.wmax)then
       SAL0=0.0d0
       SAL1=0.0d0
       SAG0=0.0d0
       SAG1=0.0d0
       SA1=0.0d0
       SPAL=0.0d0
       DO 3 I=1,3
       SPAL=SPAL+ALPHA(IQ,I,I)
       DO 3 J=1,3
       SAL0=SAL0+ALPHA(IQ,I,I)*ALPHA(IQ,J,J)
       SAL1=SAL1+ALPHA(IQ,I,J)*ALPHA(IQ,I,J)
       SAG0=SAG0+ALPHA(IQ,I,I)*    G(IQ,J,J)
       SAG1=SAG1+ALPHA(IQ,I,J)*    G(IQ,I,J)
       DO 3 K=1,3
       DO 3 L=1,3
3      SA1=SA1+EPS(I,K,L)*ALPHA(IQ,I,J)*A(IQ,K,L,J)/3.0d0
c
c      IL+IR, backscattering, DCP_I,DCP_II
       SDCP180  =6.0d0*(6.0d0*SAL1-2.0d0*SAL0)
       SDCPII180=6.0d0*(      SAL1+3.0d0*SAL0)
c
c      IL+/-IR, backscattering
       S180=7.0d0*SAL1+SAL0
c      write(6,*)ECM,SAL0,SAL1,SAG0,SAG1,SA1
       D180=8.0d0*(3.0d0*SAG1-SAG0+SA1)
       D=D180
       IF(S180.gt.1.0d-9)D=D/S180
c
c      IL+/-IR, 90o, x
       S90X=7.0d0*SAL1+SAL0
       D90X=2.0d0*(7.0d0*SAG1+SAG0+SA1)
       DX=D90X
       IF(S90X.gt.1.0d-9)DX=DX/S90X
c
c      IL+/-IR, 90o, z
       S90Z=2.0d0*(3.0d0*SAL1-SAL0)
       D90Z=4.0d0*(3.0d0*SAG1-SAG0-SA1)
       DZ=D90Z
       IF(S90Z.gt.1.0d-9)DZ=DZ/S90Z
c
c      polarized backscattering components (YDY+YDX=S180):
       YDY=3.0d0*SAL1-SAL0
       YDX=4.0d0*SAL1+2.0d0*SAL0
       if(YDX.gt.1.0d-9)then
        P1=YDY/YDX
       else
        P1=0.0d0
       endif
      
c      calculate degree of circularity
       down=        7.0d0*SAL1+ 1.0d0*SAL0
       doc=5.0d0*(  1.0d0*SAL1- 1.0d0*SAL0)
       if(down.ne.0.0d0)doc=doc/down

c      calculate depolarization
       down= 6.0d0*SAL0-12.0d0*SAL1
       ro  =-3.0d0*SAL0+ 9.0d0*SAL1
       if(down.ne.0.0d0)ro=ro/down
c
       AMU=1822.0d0
       BOHR=0.529177d0
       Alpha2=SAL0/9.0d0*(AMU*BOHR**4)
       SpA3=SPAL/3.0d0*dsqrt((AMU*BOHR**4))
       beta2=0.5d0*(3*SAL1-SAL0)*(AMU*BOHR**4)
       SDCP180  =SDCP180  *(AMU*BOHR**4)
       SDCPII180=SDCPII180*(AMU*BOHR**4)
       YDY=6.0d0*YDY*(AMU*BOHR**4)
       YDX=6.0d0*YDX*(AMU*BOHR**4)
c      gpisvejc=980.0
       gpisvejc=(AMU*BOHR**5)*1.0d4*2.0d0*4.0d0*atan(1.0d0)/EXCA
c
c      Barron's tensor invariants:
c      alpha G = [(1/3)Sp alpha] [(1/3)Sp G'] etc
c      beta(A)^2=Delta2 in Gaussian:
       alphag=SAG0/9.0d0                *gpisvejc
       betaa=SA1*3.0d0/2.0d0            *gpisvejc
       betag=(3.0d0*SAG1-SAG0)/2.0d0    *gpisvejc
       WRITE(44,4305)IQ,ECM,Alpha2,beta2,alphag,betag,betaa
4305   FORMAT(i6,f10.2,5f11.3)
c
c      DCPI 180
c      ram3=24.0d0*beta2 = 12 (3al_ab al_ab -al_aa al_bb)
       roa3=8.0d0*(12.0d0*betag+4.0d0*betaa)
c
c      polarized ROA,90deg:
c      roa90x=45.0d0*alphag+7.0d0*betag+betaa
c      ram90x=0.5d0*(45.0d0*alpha2+7.0d0*beta2)
c      depolarized ROA,90deg:
c      roa90z=6.0d0*(betag-betaa/3.0d0)
c      ram90z=0.5d0*6.0d0*beta2
c
c      unknown:
       S90X=S90X*(AMU*BOHR**4)
       D90X=D90x*gpisvejc
       DX=D90X
       IF(S90X.gt.1.0d-9)DX=DX/S90X
c
c      ICP(90)=90Z
       roa2=8.0d0*(3.0d0*betag-betaa)
       ram2=12.0d0*beta2
       CID2=0.0d0
       if(ram2.gt.1.0d-9)CID2=roa2/ram2
c
c      ICP/SCP(180):
c      8/c*(12*beta(G')^2 + 4*beta(A)^2):
       roa1=96.0d0*(betag+betaa/3.0d0)
c      ram1 = 6(7 al_ab al_ab + al_aa al_bb)
       ram1=4.0d0*(45.0d0*SpA3**2+7.0d0*beta2)
       CID1=0.0d0
       if(ram1.gt.1.0d-9)CID1=roa1/ram1

       IRR180 = 24.0d0*beta2
       IRL180 =180.0d0*SpA3**2+4.0d0*beta2


       WRITE(4,3000)IQ,ECM, YDX,YDY,ram1, CID2,DX,CID1, P1,
     1 roa1,doc,SDCP180,roa3,SDCPII180,IRR180,IRL180
3000   FORMAT(I5,f9.2,6e12.4,f6.3,e12.4,f7.3,7e12.4)
       if(lor)then
        do 6 i=1,3
        call setjk(i,j,k)
c       S0(x)=(1/2)(azz azz + ayy ayy + azy azy + azy azy)
c       S0(y)=(1/2)(axx axx + azz azz + axz axz + axz axz)
c       S0(z)=(1/2)(axx axx + ayy ayy + axy axy + axy axy)
c                    jj  jj    kk  kk    jk  jk    jk  jk
        ram1=180.0d0*0.5d0*AMU*BOHR**4*(
     1  ALPHA(IQ,j,j)*ALPHA(IQ,j,j)+ALPHA(IQ,k,k)*ALPHA(IQ,k,k)+
     1  ALPHA(IQ,j,k)*ALPHA(IQ,j,k)+ALPHA(IQ,j,k)*ALPHA(IQ,j,k))
c      -2ayx(Gxy+Gyx)+(ayy-axx)(Gxx-Gyy)
c         kj  jk  kj    kk  jj   jj  kk
c       +(2w/3)(ayyAx zy-axxAy zx+ayxAx zx-axyAy zy)
c                kk j ik  jj k ij  kj j ij  jk k ik
        roa1=-180.0d0*gpisvejc*(
     1  -2.0d0*ALPHA(IQ,j,k)*(G(IQ,j,k)+G(IQ,k,j))
     1  +(ALPHA(IQ,k,k)-ALPHA(IQ,j,j))*(G(IQ,j,j)-G(IQ,k,k))
     1  +2.0d0/3.0d0*
     1  (ALPHA(IQ,k,k)*A(IQ,j,i,k)-ALPHA(IQ,j,j)*A(IQ,k,i,j)
     1  +ALPHA(IQ,k,j)*A(IQ,j,i,j)-ALPHA(IQ,j,k)*A(IQ,k,i,k)))
6       WRITE(9+i,3000)IQ,ECM,YDX,YDY,ram1,CID2,DX,CID1,P1,roa1
       endif
c
      endif
      endif
      ENDIF

1     CONTINUE

      call closefiles(lor)

      RETURN
      END

      SUBROUTINE RSi(N,ALr,ALi,Gr,Gi,Ar,Ai,EXCA,wmin,wmax,
     1E,luseA,luseG,lor,IQ)
c
c     supposedly G/W is supplied
c
      IMPLICIT none
      integer*4 IQ,I,J,K,L,N
      real*8 ALr(N,3,3),Gr(N,3,3),Ar(N,3,3,3),E(*),EXCA,wmin,wmax,
     1ALi(N,3,3),Gi(N,3,3),Ai(N,3,3,3),
     1CM,ECM,SAL0,SAL1,SAG0,SAG1,SA1,SDCP180,SDCPII180,
     1S180,D180,D,S90X,D90X,DX,S90Z,D90Z,DZ,YDY,YDX,P1,down,doc,
     1ro,AMU,BOHR,beta2,gpisvejc,alphag,betaa,
     1betag,alpha2,roa3,roa2,ram2,CID2,roa1,
     1ram1,EPS,CID1
      logical luseA,luseG,lor
c
      call init(N,CM,luseA,luseG,Ar,Gr)
      call init(N,CM,luseA,luseG,Ai,Gi)
c
      ECM=E(IQ)*CM
      IF(ABS(ECM).GT.0.01d0)THEN
      if(wmin.eq.0.0d0.or.ECM.gt.wmin)then
      if(wmax.eq.0.0d0.or.ECM.lt.wmax)then
       SAL0=0.0d0
       SAL1=0.0d0
       SAG0=0.0d0
       SAG1=0.0d0
       SA1=0.0d0
       DO 3 I=1,3
       DO 3 J=1,3
       SAL0=SAL0+ALr(IQ,I,I)*ALr(IQ,J,J)+ALi(IQ,I,I)*ALi(IQ,J,J)
       SAL1=SAL1+ALr(IQ,I,J)*ALr(IQ,I,J)+ALi(IQ,I,J)*ALi(IQ,I,J)

c      - Im(alpha G*)
       SAG0=SAG0+ALr(IQ,I,I)* Gi(IQ,J,J)-ALi(IQ,I,I)* Gr(IQ,J,J)
       SAG1=SAG1+ALr(IQ,I,J)* Gi(IQ,I,J)-ALi(IQ,I,J)* Gr(IQ,I,J)

       DO 3 K=1,3
       DO 3 L=1,3
3      SA1=SA1+EPS(I,K,L)*( ALr(IQ,I,J)*Ar(IQ,K,L,J)
     1                     +ALi(IQ,I,J)*Ai(IQ,K,L,J))/3.0d0
c
c      IL+IR, backscattering, DCP_I,DCP_II
       SDCP180  =6.0d0*(6.0d0*SAL1-2.0d0*SAL0)
       SDCPII180=6.0d0*(      SAL1+3.0d0*SAL0)
c
c      IL+/-IR, backscattering
       S180=7.0d0*SAL1+SAL0
c      write(6,*)ECM,SAL0,SAL1,SAG0,SAG1,SA1
       D180=8.0d0*(3.0d0*SAG1-SAG0+SA1)
       D=D180
       IF(S180.gt.1.0d-9)D=D/S180
c
c      IL+/-IR, 90o, x
       S90X=7.0d0*SAL1+SAL0
       D90X=2.0d0*(7.0d0*SAG1+SAG0+SA1)
       DX=D90X
       IF(S90X.gt.1.0d-9)DX=DX/S90X
c
c      IL+/-IR, 90o, z
       S90Z=2.0d0*(3.0d0*SAL1-SAL0)
       D90Z=4.0d0*(3.0d0*SAG1-SAG0-SA1)
       DZ=D90Z
       IF(S90Z.gt.1.0d-9)DZ=DZ/S90Z
c
c      polarized backscattering components (YDY+YDX=S180):
       YDY=3.0d0*SAL1-SAL0
       YDX=4.0d0*SAL1+2.0d0*SAL0
       if(YDX.gt.1.0d-9)then
        P1=YDY/YDX
       else
        P1=0.0d0
       endif
      
c      calculate degree of circularity
       down=        7.0d0*SAL1+ 1.0d0*SAL0
       doc=5.0d0*(  1.0d0*SAL1- 1.0d0*SAL0)
       if(down.ne.0.0d0)doc=doc/down

c      calculate depolarization
       down= 6.0d0*SAL0-12.0d0*SAL1
       ro  =-3.0d0*SAL0+ 9.0d0*SAL1
       if(down.ne.0.0d0)ro=ro/down
c
       AMU=1822.0d0
       BOHR=0.529177d0
       alpha2=SAL0/9.0d0*(AMU*BOHR**4)
       beta2=0.5d0*(3*SAL1-SAL0)*(AMU*BOHR**4)
       SDCP180  =SDCP180  *(AMU*BOHR**4)
       SDCPII180=SDCPII180*(AMU*BOHR**4)
       YDY=6.0d0*YDY*(AMU*BOHR**4)
       YDX=6.0d0*YDX*(AMU*BOHR**4)
c      gpisvejc=980.0
       gpisvejc=(AMU*BOHR**5)*1.0d4*2.0d0*4.0d0*atan(1.0d0)/EXCA
c
c      Barron's tensor invariants:
c      alpha G = [(1/3)Sp alpha] [(1/3)Sp G'] etc
c      beta(A)^2=Delta2 in Gaussian:
       alphag=SAG0/9.0d0                *gpisvejc
       betag=(3.0d0*SAG1-SAG0)/2.0d0    *gpisvejc
       betaa=SA1*3.0d0/2.0d0            *gpisvejc

       WRITE(44,4305)IQ,ECM,alpha2,beta2,alphag,betag,betaa
4305   FORMAT(i6,f10.2,5f11.3)
c
c      DCPI 180
c      ram3=24.0d0*beta2 = 12 (3al_ab al_ab -al_aa al_bb)
       roa3=8.0d0*(12.0d0*betag+4.0d0*betaa)
c
c      polarized ROA,90deg:
c      roa90x=45.0d0*alphag+7.0d0*betag+betaa
c      ram90x=0.5d0*(45.0d0*alpha2+7.0d0*beta2)
c      depolarized ROA,90deg:
c      roa90z=6.0d0*(betag-betaa/3.0d0)
c      ram90z=0.5d0*6.0d0*beta2
c
c      unknown:
       S90X=S90X*(AMU*BOHR**4)
       D90X=D90x*gpisvejc
       DX=D90X
       IF(S90X.gt.1.0d-9)DX=DX/S90X
c
c      ICP(90)=90Z
       roa2=8.0d0*(3.0d0*betag-betaa)
       ram2=12.0d0*beta2
       CID2=0.0d0
       if(ram2.gt.1.0d-9)CID2=roa2/ram2
c
c      ICP/SCP(180):
c      8/c*(12*beta(G')^2 + 4*beta(A)^2):
       roa1=96.0d0*(betag+betaa/3.0d0)
c      ram1 = 6(7 al_ab al_ab + al_aa al_bb)
       ram1=4.0d0*(45.0d0*alpha2+7.0d0*beta2)
       CID1=0.0d0
       if(ram1.gt.1.0d-9)CID1=roa1/ram1

       WRITE(4,3000)IQ,ECM, YDX,YDY,ram1, CID2,DX,CID1, P1,
     1 roa1,doc,SDCP180,roa3,SDCPII180
3000   FORMAT(I5,f9.2,6e12.4,f6.3,e12.4,f7.3,5e12.4)
       if(lor)then
        do 6 i=1,3
        call setjk(i,j,k)
c       S0(x)=(1/2)(azz azz + ayy ayy + azy azy + azy azy)
c       S0(y)=(1/2)(axx axx + azz azz + axz axz + axz axz)
c       S0(z)=(1/2)(axx axx + ayy ayy + axy axy + axy axy)
c                    jj  jj    kk  kk    jk  jk    jk  jk
        ram1=180.0d0*0.5d0*AMU*BOHR**4*(
     1  ALr(IQ,j,j)*ALr(IQ,j,j)+ALr(IQ,k,k)*ALr(IQ,k,k)+
     1  ALr(IQ,j,k)*ALr(IQ,j,k)+ALr(IQ,j,k)*ALr(IQ,j,k)+
     1  ALi(IQ,j,j)*ALi(IQ,j,j)+ALi(IQ,k,k)*ALi(IQ,k,k)+
     1  ALi(IQ,j,k)*ALi(IQ,j,k)+ALi(IQ,j,k)*ALi(IQ,j,k))
c      -2ayx(Gxy+Gyx)+(ayy-axx)(Gxx-Gyy)
c         kj  jk  kj    kk  jj   jj  kk
c       +(2w/3)(ayyAx zy-axxAy zx+ayxAx zx-axyAy zy)
c                kk j ik  jj k ij  kj j ij  jk k ik
        roa1=-180.0d0*gpisvejc*(
     1  -2.0d0*ALr(IQ,j,k)*(Gr(IQ,j,k)+Gr(IQ,k,j))
     1  +(ALr(IQ,k,k)-ALr(IQ,j,j))*(Gr(IQ,j,j)-Gr(IQ,k,k))
     1  -2.0d0*ALi(IQ,j,k)*(Gi(IQ,j,k)+Gi(IQ,k,j))
     1  +(ALi(IQ,k,k)-ALi(IQ,j,j))*(Gi(IQ,j,j)-Gi(IQ,k,k))
     1  +2.0d0/3.0d0*
     1  (ALr(IQ,k,k)*Ar(IQ,j,i,k)-ALr(IQ,j,j)*Ar(IQ,k,i,j)
     1  +ALr(IQ,k,j)*Ar(IQ,j,i,j)-ALr(IQ,j,k)*Ar(IQ,k,i,k))
     1  +2.0d0/3.0d0*
     1  (ALi(IQ,k,k)*Ai(IQ,j,i,k)-ALi(IQ,j,j)*Ai(IQ,k,i,j)
     1  +ALi(IQ,k,j)*Ai(IQ,j,i,j)-ALi(IQ,j,k)*Ai(IQ,k,i,k)))
6       WRITE(9+i,3000)IQ,ECM,YDX,YDY,ram1,CID2,DX,CID1,P1,roa1
       endif
c
      endif
      endif
      ENDIF

      RETURN
      END

      SUBROUTINE TRTEN(S,N,M,ALPHA,A,G)
      IMPLICIT REAL*8 (A-H,O-Z)
      implicit integer*4 (i-n)
      DIMENSION S(N,N),A(N,3,3,3),ALPHA(N,3,3),G(N,3,3)
      real*8,allocatable::V(:)
      allocate(V(N))
      DO 1 I=1,3
      DO 1 J=1,3
      DO 3 JQ=1,N
3     V(JQ)=ALPHA(JQ,I,J)
      DO 4 IQ=1,M
      ALPHA(IQ,I,J)=0.0d0
      DO 4 IX=1,N
4     ALPHA(IQ,I,J)=ALPHA(IQ,I,J)+S(IX,IQ)*V(IX)
      DO 5 JQ=1,N
5     V(JQ)=G(JQ,I,J)
      DO 6 IQ=1,M
      G(IQ,I,J)=0.0d0
      DO 6 IX=1,N
6     G(IQ,I,J)=G(IQ,I,J)+S(IX,IQ)*V(IX)
      DO 1 K=1,3
      DO 7 JQ=1,N
7     V(JQ)=A(JQ,I,J,K)
      DO 1 IQ=1,M
      A(IQ,I,J,K)=0.0D0
      DO 1 IX=1,N
1     A(IQ,I,J,K)=A(IQ,I,J,K)+S(IX,IQ)*V(IX)
      RETURN
      END

      SUBROUTINE DORAMANi(N,M,ALr,ALi,Gr,Gi,Ar,Ai,EXCA,wmin,wmax,
     1E,luseA,luseG,lor)
c
c     supposedly G/W is supplied
c
      IMPLICIT none
      integer*4 IQ,I,J,K,L,N,M
      real*8 ALr(N,3,3),Gr(N,3,3),Ar(N,3,3,3),E(*),EXCA,wmin,wmax,
     1ALi(N,3,3),Gi(N,3,3),Ai(N,3,3,3),
     1CM,ECM,SAL0,SAL1,SAG0,SAG1,SA1,SDCP180,SDCPII180,
     1S180,D180,D,S90X,D90X,DX,S90Z,D90Z,DZ,YDY,YDX,P1,down,doc,
     1ro,AMU,BOHR,beta2,gpisvejc,alphag,betaa,
     1betag,alpha2,roa3,roa2,ram2,CID2,roa1,
     1ram1,EPS,CID1
      logical luseA,luseG,lor
c
      call init(N,CM,luseA,luseG,Ar,Gr)
      call init(N,CM,luseA,luseG,Ai,Gi)
      call openfiles(lor)
c
      DO 1 IQ=1,M
      ECM=E(IQ)*CM
      IF(ABS(ECM).GT.0.01d0)THEN
      if(wmin.eq.0.0d0.or.ECM.gt.wmin)then
      if(wmax.eq.0.0d0.or.ECM.lt.wmax)then
       SAL0=0.0d0
       SAL1=0.0d0
       SAG0=0.0d0
       SAG1=0.0d0
       SA1=0.0d0
       DO 3 I=1,3
       DO 3 J=1,3
       SAL0=SAL0+ALr(IQ,I,I)*ALr(IQ,J,J)+ALi(IQ,I,I)*ALi(IQ,J,J)
       SAL1=SAL1+ALr(IQ,I,J)*ALr(IQ,I,J)+ALi(IQ,I,J)*ALi(IQ,I,J)

c      - Im(alpha G*)
       SAG0=SAG0+ALr(IQ,I,I)* Gi(IQ,J,J)-ALi(IQ,I,I)* Gr(IQ,J,J)
       SAG1=SAG1+ALr(IQ,I,J)* Gi(IQ,I,J)-ALi(IQ,I,J)* Gr(IQ,I,J)

       DO 3 K=1,3
       DO 3 L=1,3
3      SA1=SA1+EPS(I,K,L)*( ALr(IQ,I,J)*Ar(IQ,K,L,J)
     1                     +ALi(IQ,I,J)*Ai(IQ,K,L,J))/3.0d0
c
c      IL+IR, backscattering, DCP_I,DCP_II
       SDCP180  =6.0d0*(6.0d0*SAL1-2.0d0*SAL0)
       SDCPII180=6.0d0*(      SAL1+3.0d0*SAL0)
c
c      IL+/-IR, backscattering
       S180=7.0d0*SAL1+SAL0
c      write(6,*)ECM,SAL0,SAL1,SAG0,SAG1,SA1
       D180=8.0d0*(3.0d0*SAG1-SAG0+SA1)
       D=D180
       IF(S180.gt.1.0d-9)D=D/S180
c
c      IL+/-IR, 90o, x
       S90X=7.0d0*SAL1+SAL0
       D90X=2.0d0*(7.0d0*SAG1+SAG0+SA1)
       DX=D90X
       IF(S90X.gt.1.0d-9)DX=DX/S90X
c
c      IL+/-IR, 90o, z
       S90Z=2.0d0*(3.0d0*SAL1-SAL0)
       D90Z=4.0d0*(3.0d0*SAG1-SAG0-SA1)
       DZ=D90Z
       IF(S90Z.gt.1.0d-9)DZ=DZ/S90Z
c
c      polarized backscattering components (YDY+YDX=S180):
       YDY=3.0d0*SAL1-SAL0
       YDX=4.0d0*SAL1+2.0d0*SAL0
       if(YDX.gt.1.0d-9)then
        P1=YDY/YDX
       else
        P1=0.0d0
       endif
      
c      calculate degree of circularity
       down=        7.0d0*SAL1+ 1.0d0*SAL0
       doc=5.0d0*(  1.0d0*SAL1- 1.0d0*SAL0)
       if(down.ne.0.0d0)doc=doc/down

c      calculate depolarization
       down= 6.0d0*SAL0-12.0d0*SAL1
       ro  =-3.0d0*SAL0+ 9.0d0*SAL1
       if(down.ne.0.0d0)ro=ro/down
c
       AMU=1822.0d0
       BOHR=0.529177d0
       alpha2=SAL0/9.0d0*(AMU*BOHR**4)
       beta2=0.5d0*(3*SAL1-SAL0)*(AMU*BOHR**4)
       SDCP180  =SDCP180  *(AMU*BOHR**4)
       SDCPII180=SDCPII180*(AMU*BOHR**4)
       YDY=6.0d0*YDY*(AMU*BOHR**4)
       YDX=6.0d0*YDX*(AMU*BOHR**4)
c      gpisvejc=980.0
       gpisvejc=(AMU*BOHR**5)*1.0d4*2.0d0*4.0d0*atan(1.0d0)/EXCA
c
c      Barron's tensor invariants:
c      alpha G = [(1/3)Sp alpha] [(1/3)Sp G'] etc
c      beta(A)^2=Delta2 in Gaussian:
       alphag=SAG0/9.0d0                *gpisvejc
       betag=(3.0d0*SAG1-SAG0)/2.0d0    *gpisvejc
       betaa=SA1*3.0d0/2.0d0            *gpisvejc

       WRITE(44,4305)IQ,ECM,alpha2,beta2,alphag,betag,betaa
4305   FORMAT(i6,f10.2,5f11.3)
c
c      DCPI 180
c      ram3=24.0d0*beta2 = 12 (3al_ab al_ab -al_aa al_bb)
       roa3=8.0d0*(12.0d0*betag+4.0d0*betaa)
c
c      polarized ROA,90deg:
c      roa90x=45.0d0*alphag+7.0d0*betag+betaa
c      ram90x=0.5d0*(45.0d0*alpha2+7.0d0*beta2)
c      depolarized ROA,90deg:
c      roa90z=6.0d0*(betag-betaa/3.0d0)
c      ram90z=0.5d0*6.0d0*beta2
c
c      unknown:
       S90X=S90X*(AMU*BOHR**4)
       D90X=D90x*gpisvejc
       DX=D90X
       IF(S90X.gt.1.0d-9)DX=DX/S90X
c
c      ICP(90)=90Z
       roa2=8.0d0*(3.0d0*betag-betaa)
       ram2=12.0d0*beta2
       CID2=0.0d0
       if(ram2.gt.1.0d-9)CID2=roa2/ram2
c
c      ICP/SCP(180):
c      8/c*(12*beta(G')^2 + 4*beta(A)^2):
       roa1=96.0d0*(betag+betaa/3.0d0)
c      ram1 = 6(7 al_ab al_ab + al_aa al_bb)
       ram1=4.0d0*(45.0d0*alpha2+7.0d0*beta2)
       CID1=0.0d0
       if(ram1.gt.1.0d-9)CID1=roa1/ram1

       WRITE(4,3000)IQ,ECM, YDX,YDY,ram1, CID2,DX,CID1, P1,
     1 roa1,doc,SDCP180,roa3,SDCPII180
3000   FORMAT(I5,f9.2,6e12.4,f6.3,e12.4,f7.3,5e12.4)
       if(lor)then
        do 6 i=1,3
        call setjk(i,j,k)
c       S0(x)=(1/2)(azz azz + ayy ayy + azy azy + azy azy)
c       S0(y)=(1/2)(axx axx + azz azz + axz axz + axz axz)
c       S0(z)=(1/2)(axx axx + ayy ayy + axy axy + axy axy)
c                    jj  jj    kk  kk    jk  jk    jk  jk
        ram1=180.0d0*0.5d0*AMU*BOHR**4*(
     1  ALr(IQ,j,j)*ALr(IQ,j,j)+ALr(IQ,k,k)*ALr(IQ,k,k)+
     1  ALr(IQ,j,k)*ALr(IQ,j,k)+ALr(IQ,j,k)*ALr(IQ,j,k)+
     1  ALi(IQ,j,j)*ALi(IQ,j,j)+ALi(IQ,k,k)*ALi(IQ,k,k)+
     1  ALi(IQ,j,k)*ALi(IQ,j,k)+ALi(IQ,j,k)*ALi(IQ,j,k))
c      -2ayx(Gxy+Gyx)+(ayy-axx)(Gxx-Gyy)
c         kj  jk  kj    kk  jj   jj  kk
c       +(2w/3)(ayyAx zy-axxAy zx+ayxAx zx-axyAy zy)
c                kk j ik  jj k ij  kj j ij  jk k ik
        roa1=-180.0d0*gpisvejc*(
     1  -2.0d0*ALr(IQ,j,k)*(Gr(IQ,j,k)+Gr(IQ,k,j))
     1  +(ALr(IQ,k,k)-ALr(IQ,j,j))*(Gr(IQ,j,j)-Gr(IQ,k,k))
     1  -2.0d0*ALi(IQ,j,k)*(Gi(IQ,j,k)+Gi(IQ,k,j))
     1  +(ALi(IQ,k,k)-ALi(IQ,j,j))*(Gi(IQ,j,j)-Gi(IQ,k,k))
     1  +2.0d0/3.0d0*
     1  (ALr(IQ,k,k)*Ar(IQ,j,i,k)-ALr(IQ,j,j)*Ar(IQ,k,i,j)
     1  +ALr(IQ,k,j)*Ar(IQ,j,i,j)-ALr(IQ,j,k)*Ar(IQ,k,i,k))
     1  +2.0d0/3.0d0*
     1  (ALi(IQ,k,k)*Ai(IQ,j,i,k)-ALi(IQ,j,j)*Ai(IQ,k,i,j)
     1  +ALi(IQ,k,j)*Ai(IQ,j,i,j)-ALi(IQ,j,k)*Ai(IQ,k,i,k)))
6       WRITE(9+i,3000)IQ,ECM,YDX,YDY,ram1,CID2,DX,CID1,P1,roa1
       endif
c
      endif
      endif
      ENDIF
1     CONTINUE
      call closefiles(lor)
      RETURN
      END
c     ====================
      subroutine readpol(a)
      real*8 a(3,3)
      open(90,file='POL.TTT')
      read(90,*)
      read(90,*)
      read(90,*)a(1,1),a(1,2),a(2,2),a(1,3),a(2,3),a(3,3)
      close(90)
      a(2,1)=a(1,2)
      a(3,1)=a(1,3)
      a(3,2)=a(2,3)
      return
      end
c     ====================
      SUBROUTINE rules(ALPHA,A,G,NAT)
C     project translations/rotations
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION ALPHA(3*NAT,3,3),G(3*NAT,3,3),A(3*NAT,3,3,3),
     1ap(3,3),aap(3,3,3),APOL(3,3)
      character*2 xyz(3)
      logical lpol
      data xyz/' X',' Y',' Z'/
c
      write(6,6077)' Sum rules-alpha:'
6077  format(a17)
      do 11 ii=1,3
      do 21 ix=1,3
      do 21 jx=1,3
      ap(ix,jx)=0.0d0
      do 21 ia=2,NAT
21    ap(ix,jx)=ap(ix,jx)+ALPHA((ia-1)*3+ii,ix,jx)
      write(6,6017)xyz(ii)
6017  format(a2)
      write(6,600)((  -ap(ix,jx),jx=1,3),ix=1,3)
600   format(9F8.2)
11    write(6,600)((ALPHA(ii,ix,jx),jx=1,3),ix=1,3)

      inquire(file='POL.TTT',exist=lpol)
      if(lpol)then
       call readpol(APOL)
       write(6,*)'Static polarizability read in'
      else
       do 61 ix=1,3
       do 61 jx=1,3
61     APOL(ix,jx)=0.0d0
      endif

      write(6,6077)' Sum rules-G    :'
      do 1 ii=1,3
      do 2 ix=1,3
      do 2 jx=1,3
      ap(ix,jx)=0.0d0
      do 2 ia=2,NAT
2     ap(ix,jx)=ap(ix,jx)+    G((ia-1)*3+ii,ix,jx)
      write(6,6017)xyz(ii)
      write(6,600)((-ap(ix,jx),jx=1,3),ix=1,3)
      do 3 ix=1,3
      do 3 jx=1,3
      ap(ix,jx)=0.0d0
      do 3 kx=1,3
3     ap(ix,jx)=ap(ix,jx)+0.5d0*EPS(jx,ii,kx)*APOL(ix,kx)
1     write(6,600)((G(ii,ix,jx)+ap(ix,jx),jx=1,3),ix=1,3)
c
      write(6,6077)' Sum rules-A    :'
      do 14 ii=1,3
      do 12 ix=1,3
      do 12 jx=1,3
      do 12 kx=1,3
      aap(ix,jx,kx)=0.0d0
      do 12 ia=2,NAT
12    aap(ix,jx,kx)=aap(ix,jx,kx)+A(3*(ia-1)+ii,ix,jx,kx)
      do 14 kx=1,3
      write(6,6018)xyz(ii),xyz(kx)
6018  format(2a2)
      write(6,600)((-aap(ix,jx,kx),jx=1,3),ix=1,3)
      do 41 ix=1,3
      do 41 jx=1,3
      ap(ix,jx)=0.0d0
      if(jx.eq.kx)ap(ix,jx)=ap(ix,jx)      +APOL(ix,ii)
      if(ii.eq.jx)ap(ix,jx)=ap(ix,jx)-1.5d0*APOL(ix,kx)
41    if(ii.eq.kx)ap(ix,jx)=ap(ix,jx)-1.5d0*APOL(ix,jx)
14    write(6,600)((A(ii,ix,jx,kx)+ap(ix,jx),jx=1,3),ix=1,3)

      return
      end
c     ====================
      SUBROUTINE TTTT(N,ALPHA,A,G,NAT,ICO,R)
c
C     ICO = 1: Transforms A and G to local origins
C     ICO = 2: Transforms A and G to common origin
C     ICO = 3: Sets A and G to zero 
c     ICO = 4: Delete G
c     ICO = 5: Delete A
C     R supposed to be passed in AU
C
      IMPLICIT none
      INTEGER*4 N,NAT,ICO,I,J,L,IE,II,K,M,IA,IB,IC,ID
      real*8 ALPHA(N,3,3),G(N,3,3),A(N,3,3,3),R(3,N/2),
     1S,EPS,SUM
      real*8 ,allocatable::x(:,:)
      allocate(X(3,N/3))
C
      DO 6 I=1,3
      DO 6 J=1,NAT
6     X(I,J)=R(I,J)
C
      IF(ICO.EQ.1)then
       S= 1.0d0
       WRITE(6,*)' A,G transformed into atomic origins'
      elseif(ICO.EQ.2)then
       S=-1.0d0
       WRITE(6,*)' A,G transformed into common origin'
      elseif(ICO.ge.3.and.ICO.LE.5)then
       IF(ICO.eq.3.or.ICO.EQ.4)then
        do 5 II=1,N
        do 5 I=1,3
        do 5 J=1,3
5       G(II,I,J)=0.0d0
        WRITE(6,*)' G set to zero'
       endif
       IF(ICO.eq.3.or.ICO.EQ.5)then
        do 7 II=1,N
        do 7 I=1,3
        do 7 J=1,3
        do 7 K=1,3
7       A(II,I,J,K)=0.0d0
        WRITE(6,*)' A set to zero'
       endif
       return
      else
       write(6,*)ICO
       call report('TTTT unknown instruction')
      endif
C
      DO 1 I=1,3
      DO 1 J=1,3
      DO 1 L=1,NAT
      DO 1 IE=1,3
      II=3*(L-1)+IE
      DO 1 K=1,3
      DO 1 M=1,3
c     G-second index-magnetic
1     G(II,I,J)=G(II,I,J)+S*0.5d0*EPS(J,K,M)*X(K,L)*ALPHA(II,I,M)
C
      DO 3 IA=1,3
      DO 3 IB=1,3
      DO 3 IC=1,3
      DO 3 L=1,NAT
      DO 3 IE=1,3
      II=3*(L-1)+IE
      SUM=0.0d0
      IF(IB.EQ.IC)THEN
       DO 4 ID=1,3
4      SUM=SUM+X(ID,L)*ALPHA(II,IA,ID)
      ENDIF
      A(II,IA,IB,IC)=A(II,IA,IB,IC)+SUM*S
     1-1.5d0*S*(X(IB,L)*ALPHA(II,IA,IC)+X(IC,L)*ALPHA(II,IA,IB))
3     CONTINUE
c
      RETURN
      END
c     ====================================
      subroutine setjk(i,j,k)
      integer*4 i,j,k
      j=i+1
      if(j.gt.3)j=1
      k=j+1
      if(k.gt.3)k=1
      return
      end

      function inw(s)
      integer*4 nst,inw
      character*80 filename
      character*(*)s
      nst=0
      open(9,file=s,status='old')
1     read(9,90,end=99,err=99)filename
90    format(a80)
c     avoid compiler messages:
      if(filename.eq.' ')filename=' '
      nst=nst+1
      goto 1
   99 close(9)
      write(6,601)nst
 601  format(i6,' lines in ',$)
      write(6,*)s
      inw=nst
      return
      end

      subroutine fille(we,nw)
      implicit none
      integer*4 nw,i
      real*8 we(nw)
      character*80 fn,s
      open(8,file='TTT.LST')
      do 1 i=1,nw
      read(8,80)fn
 80   format(a80)
      open(9,file=fn,status='old')
      read(9,*)
      read(9,80)s
      read(s(21:31),*)we(i)
1     close(9)
      close(8)
      return
      end

      subroutine report(s)
      character*(*) s
      write(6,*)s
      stop
      end

      subroutine mkt(N,nw,wd,we,P,A,G,Pi,Ai,Gi,lcompl,iq,e,na,iwr)
      implicit none
      integer*4 nw,i,j,k,l,i1,i2,N,NAT,iq,na,m,q,iwr
      real*8 wd,we(nw),w1,w2,d,P(N,3,3),A(N,3,3,3),G(N,3,3),
     1Pi(N,3,3),Ai(N,3,3,3),Gi(N,3,3),CM,e
      real*8,allocatable::P1(:,:,:),G1(:,:,:),A1(:,:,:,:),P2(:,:,:),
     1G2(:,:,:),A2(:,:,:,:),P1i(:,:,:),G1i(:,:,:),A1i(:,:,:,:),
     1P2i(:,:,:),G2i(:,:,:),A2i(:,:,:,:),am(:,:),vm(:),ami(:,:),
     1PX(:,:,:,:),GX(:,:,:,:),AX(:,:,:,:,:),
     1PXi(:,:,:,:),GXi(:,:,:,:),AXi(:,:,:,:,:)
      logical lcompl
      character*80 fn
      allocate(P1(N,3,3),A1(N,3,3,3),G1(N,3,3),
     1P1i(N,3,3),A1i(N,3,3,3),G1i(N,3,3),
     1P2 (N,3,3),A2 (N,3,3,3),G2 (N,3,3),
     1P2i(N,3,3),A2i(N,3,3,3),G2i(N,3,3))
      if(na.gt.0)then
       m=na+1
       if(nw.lt.m)call report('nw < na+1')
       allocate(am(m,m),vm(m),ami(m,m))
       am=0.0d0
       do 4 i=1,m
       do 4 j=1,m
       do 4 k=1,nw
4      am(i,j)=am(i,j)+we(k)**(i+j-2)
       call inv(am,ami,m)
       allocate(PX(nw,N,3,3),AX(nw,N,3,3,3),GX(nw,N,3,3),
     1 PXi(nw,N,3,3),AXi(nw,N,3,3,3),GXi(nw,N,3,3))
       open(8,file='TTT.LST')
       do 5 i=1,nw
       read(8,80)fn
       call READTTT(N,NAT,P1,G1,A1,0,fn)
       if(lcompl)then
        call READTTT(N,NAT,P1i,G1i,A1i,3,fn)
        do 6 l=1,N
        do 6 j=1,3
        do 6 k=1,3
        PXi(i,l,j,k)=P1i(l,j,k)
        GXi(i,l,j,k)=G1i(l,j,k)
        do 6 q=1,3
 6      AXi(i,l,j,k,q)=A1i(l,j,k,q)
       endif
       do 5 l=1,N
       do 5 j=1,3
       do 5 k=1,3
       PX(i,l,j,k)=P1(l,j,k)
       GX(i,l,j,k)=G1(l,j,k)
       do 5 q=1,3
 5     AX(i,l,j,k,q)=A1(l,j,k,q)
       close(8)

       call mxx(N,m,nw,vm,PX,P,ami,we,wd,iwr)
       call mxx(N,m,nw,vm,GX,G,ami,we,wd,iwr)
       call mx3(N,m,nw,vm,AX,A,ami,we,wd,iwr)
       if(lcompl)then
        call mxx(N,m,nw,vm,PXi,Pi,ami,we,wd,iwr)
        call mxx(N,m,nw,vm,GXi,Gi,ami,we,wd,iwr)
        call mx3(N,m,nw,vm,AXi,Ai,ami,we,wd,iwr)
       endif

      else
       if(nw.lt.2)call report('not enough frequencies')

       NAT=N/3
       CM=219470.0d0
c      closest frequency:
       w1=we(1)
       i1=1
       d=dabs(we(1)-wd)
       do 1 i=2,nw
       if(dabs(we(i)-wd).lt.d)then
        d=dabs(we(i)-wd)
        w1=we(i)
        i1=i
       endif
1      continue
c      second closest frequency:
       i2=1
       if(i1.eq.1)i2=2
       w2=we(i2)
       d=dabs(w2-wd)
       do 2 i=2,nw
       if(dabs(we(i)-wd).lt.d.and.i.ne.i1)then
        d=dabs(we(i)-wd)
        w2=we(i)
        i2=i
       endif
2      continue
       P1i=0.0d0
       P2i=0.0d0
       G1i=0.0d0
       G2i=0.0d0
       A1i=0.0d0
       A2i=0.0d0
c      load the two closest derivatives
       open(8,file='TTT.LST')
       do 3 i=1,nw
       read(8,80)fn
80     format(A80)
       if(i.eq.i1)then
        call READTTT(N,NAT,P1,G1,A1,0,fn)
        if(lcompl)call READTTT(N,NAT,P1i,G1i,A1i,3,fn)
       endif
       if(i.eq.i2)then
        call READTTT(N,NAT,P2,G2,A2,0,fn)
        if(lcompl)call READTTT(N,NAT,P2i,G2i,A2i,3,fn)
       endif
3      continue
       close(8)
c      approximate for wd:
       P=P1+(P2-P1)*(wd-w1)/(w2-w1)
       G=G1+(G2-G1)*(wd-w1)/(w2-w1)
       A=A1+(A2-A1)*(wd-w1)/(w2-w1)
       Pi=P1i+(P2i-P1i)*(wd-w1)/(w2-w1)
       Gi=G1i+(G2i-G1i)*(wd-w1)/(w2-w1)
       Ai=A1i+(A2i-A1i)*(wd-w1)/(w2-w1)
       write(6,600)iq,e*CM,1.0d7/wd/CM,1.0d7/w1/cm,1.0d7/w2/CM,i1,i2
600    format(i6,'. mode',f10.1,' cm-1; wd w1 w2 i1 i2:',3(f10.1,' nm')
     1 ,2i3)
      endif
      return
      end

      subroutine mxx(N,m,nw,vm,PX,P,ai,we,wd,iwr)
      integer*4 m,nw,N,ii,i,l,j,k,iwr
      real*8 vm(m),PX(nw,N,3,3),P(N,3,3),ai(m,m),wd,we(nw),
     1wi,wa,dw,w,pt
      real*8 ,allocatable::sa(:)
      allocate(sa(m))
      P=0.0d0
      do 7 l=1,N
      do 7 j=1,3
      do 7 k=1,3
      vm=0.0d0
      do 71 ii=1,m
      do 71 i=1,nw
71    vm(ii) =vm(ii)+PX(i,l,j,k)*we(i)**(ii-1)
      call mksa(sa,m,ai,vm)
      do 7 ii=1,m
7     P(l,j,k)=P(l,j,k)+sa(ii)*wd**(ii-1)

      if(iwr.gt.0)then
       write(6,*)' 1 1 3 component:'
       wi=wd
       wa=wd
       do 72 i=1,nw
       if(we(i).lt.wi)wi=we(i)
       if(we(i).gt.wa)wa=we(i)
72     write(6,600)i,we(i),PX(i,1,1,3)
600    format(i3,2f12.6)
       write(6,601)wd,P(1,1,3)
601    format(2f12.6)
       if(iwr.gt.1)then
        write(6,*)'Fitted curve:'
        dw=(wa-wi)/100.0d0
        w=wi-dw
        do 73 j=1,100
        w=w+dw
        vm=0.0d0
        do 74 ii=1,m
        do 74 i=1,nw
74      vm(ii) =vm(ii)+PX(i,1,1,3)*we(i)**(ii-1)
        call mksa(sa,m,ai,vm)
        pt=0.0d0
        do 75 ii=1,m
75      pt=pt+sa(ii)*w**(ii-1)
73      write(6,600)j,w,pt
       endif
      endif

      return
      end

      subroutine mx3(N,m,nw,vm,PX,P,ai,we,wd,iwr)
      integer*4 m,nw,N,ii,i,l,j,k,q,iwr
      real*8 vm(m),PX(nw,N,3,3,3),P(N,3,3,3),ai(m,m),wd,we(nw),pt,
     1wi,wa,dw,w
      real*8 ,allocatable::sa(:)
      allocate(sa(m))
      P=0.0d0
      do 7 l=1,N
      do 7 j=1,3
      do 7 k=1,3
      do 7 q=1,3
      vm=0.0d0
      do 71 ii=1,m
      do 71 i=1,nw
71    vm(ii) =vm(ii)+PX(i,l,j,k,q)*we(i)**(ii-1)
      call mksa(sa,m,ai,vm)
      do 7 ii=1,m
7     P(l,j,k,q)=P(l,j,k,q)+sa(ii)*wd**(ii-1)

      if(iwr.gt.0)then
       write(6,*)' 1 1 3 3 component:'
       wi=wd
       wa=wd
       do 72 i=1,nw
       if(we(i).lt.wi)wi=we(i)
       if(we(i).gt.wa)wa=we(i)
72     write(6,600)i,we(i),PX(i,1,1,3,3)
600    format(i3,2f12.6)
       write(6,601)wd,P(1,1,3,3)
601    format(2f12.6)
       if(iwr.gt.1)then
        write(6,*)'Fitted curve:'
        dw=(wa-wi)/100.0d0
        w=wi-dw
        do 73 j=1,100
        w=w+dw
        vm=0.0d0
        do 74 ii=1,m
        do 74 i=1,nw
74      vm(ii) =vm(ii)+PX(i,1,1,3,3)*we(i)**(ii-1)
        call mksa(sa,m,ai,vm)
        pt=0.0d0
        do 75 ii=1,m
75      pt=pt+sa(ii)*w**(ii-1)
73      write(6,600)j,w,pt
       endif
      endif
      return
      end

      subroutine mksa(sa,m,ai,vm)
      integer*4 m,ii,jj
      real*8 sa(m),ai(m,m),vm(m)
      sa=0.0d0
      do 75 ii=1,m
      do 75 jj=1,m
75    sa(ii)=sa(ii)+ai(ii,jj)*vm(jj)
      return
      end

      subroutine init(N,CM,luseA,luseG,A,G)
      IMPLICIT none
      integer*4 N
      real*8 CM
      real*8 A(N,3,3,3),G(N,3,3)
      logical luseA,luseG
      CM=219470.0d0
      if(.not.luseA)A=0.0d0
      if(.not.luseG)G=0.0d0
      return
      end

      subroutine openfiles(lor)
      character*8 OF(3)
      data OF/'ROAX.TAB','ROAY.TAB','ROAZ.TAB'/
      logical lor
      integer*4 i
      if(lor)then
       do 5 i=1,3
       OPEN(9+i,FILE=OF(i))
5      WRITE(9+i,3304)
      endif
      OPEN(4,FILE='ROA.TAB')
      WRITE(4,3304)
3304  FORMAT(
     1'MODE      FREQ    Ramx (180)  Ramy',6x,'Ramtot',11x,'CID90',
     27x,'CID-X',6x,'CID180',3x,'DEP',5x,'ROA180',5x,'DOC',
     15x,'RamDCPI',5x,'RoaDCPI',3x,'RamnDCPII',4x,'IRR(180)',
     14x,'IRL(180)',/,10x,'cm-1',/,149(1H-))
      OPEN(44,FILE='INV.TXT')
      WRITE(44,4304)
4304  FORMAT(4x,'IQ',9x,'W',5x,'Alpha2',6x,'Beta2',5x,'AlphaG',
     1'     Gamma2     Delta2')
      return
      end

      SUBROUTINE INV(a,ai,n)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      dimension ai(n,n),a(n,n)
      real*8,allocatable::e(:,:)

      TOL=1.0D-10
      allocate(e(n,2*n))
C
10000 DO 1  ii=1,n
      DO 1  jj=1,n
      e(ii,jj)=a(ii,jj)
1     e(ii,jj+n)=0.0D0
      do 13 ii=1,n
13    e(ii,ii+n)=1.0D0
c
      DO 2 ii=1,n-1
      if (ABS(e(ii,ii)).LE.TOL) then
       DO 3 io=ii+1,n
3      if (ABS(e(io,ii)).GT.TOL) goto 11
       tol=tol*0.50d0
       if(tol.gt.1.0d-20)goto 10000
       call report('cannot invert matrix')
c
11     CONTINUE
       DO 4 kk=1,2*n
       w=e(ii,kk)
       e(ii,kk)=e(io,kk)
4      e(io,kk)=w
      ENDIF
      eii=e(ii,ii)
      DO 5 jj=ii+1,n
      e1=e(jj,ii)/eii
      DO 6 kk=ii+1, 2*n
6     e(jj,kk)=e(jj,kk)-e(ii,kk)*e1
5     e(jj,ii)=0.0D0
2     CONTINUE
c
      DO 7  i2=n,2,-1
      eii=e(i2,i2)
      DO 7 j2=i2-1,1,-1
      e1=e(j2,i2)/eii
      DO 9 kk=1, n
9     e(j2,kk+n)=e(j2,kk+n)-e(i2,kk+n)*e1
7     e(j2,i2)=0.0d0
c
      DO 10 ii=1,n
      ei=1.0d0/e(ii,ii)
      DO 12  jj=1,n
12    ai(ii,jj)=e(ii,jj+n)*ei
10    CONTINUE
c
      RETURN
      END

      subroutine closefiles(lor)
      logical lor
      integer*4 i
      WRITE(4,3002)
3002  FORMAT(149(1H-))
      CLOSE(4)
      CLOSE(44)
      if(lor)then
       do 7 i=1,3
       WRITE(9+i,3002)
7      CLOSE(9+i)
      endif
      return
      end