      PROGRAM RELAX2X
     0COMPLEX A(4,4,2),B(2,2,2),EIFI,SPU(16,16),ISV(16,16),
     1W(2,2),IM,PW(2),SPQM(2,4,2,4),Z,Z01,Z02,EIEPSZI,SV(16,16),SE(16),
     2D(16,16),TEMP,SEO(16),AP(2,2,4),SP(2,2,4,2,2,4),SPQ
     0DIMENSION Q(2),TETA(2),FI(2),ETA(2),H(2),TAU(2),SSC(256),
     1SS(256),CHEM(2),ETETA(2),EPSZI(2),EFI(2),OMEGA(256)
      EQUIVALENCE (SP(1,1,1,1,1,1),SPU(1,1))
      CHARACTER*12 DATAFN,ISPFN,OSPFN
      CHARACTER*1 ANS
      CHARACTER*72 STIT
C
C     **************************    START   ***************************
C
      WRITE(*,600)
 600  FORMAT(////////////,35X,'R E L A X 2 X',/,37X,'by mm/s',/////,7X,
     1' Program for computing M”ssbauer-spectrum of relaxing iron ',
     2' nucleous',////,30X,'Press ENTER to begin',////)
      READ(*,'(A)') ANS
      ANS='N'
C
C     **************************    MENU    ***************************
C
   1  CALL MENU(ANS,DATAFN,ISPFN,OSPFN,STIT)
      IF ((ANS.EQ.'Y').OR.(ANS.EQ.'y')) GO TO 100
      OPEN(6,FILE='LPT1')
      WRITE(6,603) STIT,DATAFN,OSPFN
 603  FORMAT(//,35X,'R E L A X 2 X'///,37X,'by mm/s',/////
     1' SPECTRUM TITLE : ',A72,/,
     2' INPUT DATA-FILE : ',A12,//,
     3' OUTPUT SPECTRUM FILE : ',A12,//)
C
C     ************************** DATA INPUT ***************************
C
      OPEN(1,FILE=DATAFN,STATUS='OLD',FORM='FORMATTED')
      READ(1,620) NCHAN,REF,DELTA,H,TETA,FI,Q,ETA,ETETA,EPSZI,EFI,CHEM,
     1            TAU,REFGAM,TETAP,FIAP
 620  FORMAT(I3/,24(F12.5/),F12.5)
      CLOSE(1,STATUS='KEEP')
      WRITE(*,630) NCHAN,REF,DELTA,H,TETA,FI,Q,ETA,ETETA,EPSZI,EFI,CHEM,
     1            TAU,REFGAM,TETAP,FIAP
      WRITE(6,630) NCHAN,REF,DELTA,H,TETA,FI,Q,ETA,ETETA,EPSZI,EFI,CHEM,
     1            TAU,REFGAM,TETAP,FIAP
 630  FORMAT(/31X,17HINPUT PARAMETERS ///1X,18HNUMBER OF CHANNELS,
     1       29X,I3/1X,38HVELOCITY OF THE FIRST CHANNEL (MM/SEC),8X,
     2       G12.5/1X,22HCHANNEL WIDTH (MM/SEC),24X,G12.5///
     3       49X,7HSTATE 1,13X,7HSTATE 2//
     4       1X,20HMAGNETIC FIELD (KOE),26X,G12.5,8X,G12.5/
     5       1X,34HMAGNETIC FIELD POLAR ANGLES  THETA,12X,
     6       G12.5,8X,G12.5/30X,3HPHI,14X,G12.5,8X,G12.5//
     7       1X,41HQUADRUPOLE INTERACTION (EQVZZ/2) (MM/SEC),5X,G12.5,
     8       8X,G12.5/
     9       1X,20HASYMMETRY OF THE EFG,26X,G12.5,8X,G12.5/
     T       1X,34HEULER ANGLES OF THE EFG      THETA,12X,G12.5,
     1       8X,G12.5/30X,3HPSI,14X,G12.5,8X,G12.5/
     2       30X,3HPHI,14X,G12.5,8X,G12.5//
     3       1X,21HISOMER SHIFT (MM/SEC),25X,G12.5,8X,G12.5/
     4       1X,22HRELAXATION TIME (NSEC),24X,G12.5,8X,G12.5/
     5       1X,15HREFGAM (MM/SEC),41X,G12.5/
     6       1X,'OBSERVATION DIRECTION POLAR ANGLES ; THETA :',G12.5/
     7       1X,37X,'PHI   :',G12.5///)
C
C     *****************************************************************
C
      TAUN=141.
      TAUP=.0728254
      GAMMAFEL=0.5/TAUN/TAUP
      SQ3=SQRT(3.)
      SQ32=SQ3/2.
      IM=(0.,1.)
      ALFAP=-.00678522
      BETAP=.0118694
      DO 10 J=1,256
  10  SS(J)=0.
      AP(1,2,1)=(0.,0.)
      AP(1,1,4)=(0.,0.)
      AP(2,2,1)=(0.,0.)
      AP(2,1,4)=(0.,0.)
      AP(1,1,1)=CEXP(CMPLX(0.,FIAP))*0.5*(1+COS(TETAP))
      AP(2,1,1)=CEXP(CMPLX(0.,FIAP))*0.5*(1-COS(TETAP))
      AP(1,2,4)=CEXP(CMPLX(0.,-FIAP))*0.5*(1-COS(TETAP))
      AP(2,2,4)=CEXP(CMPLX(0.,-FIAP))*0.5*(1+COS(TETAP))
      AP(1,2,2)=AP(1,1,1)/SQ3
      AP(2,2,2)=AP(2,1,1)/SQ3
      AP(1,1,3)=AP(1,2,4)/SQ3
      AP(2,1,3)=AP(2,2,4)/SQ3
      AP(2,2,3)=CMPLX(SIN(TETAP)/SQ3,0.)
      AP(2,1,2)=AP(2,2,3)
      AP(1,1,2)=-AP(2,2,3)
      AP(1,2,3)=-AP(2,2,3)
C
C     ************** OBTAIN HAMILTONIAN MATRICES (A AND B) ************
C
      DO 20 I=1,2
      AST=SIN(TETA(I))
      ACT=COS(TETA(I))
      EIFI=CEXP(CMPLX(0.,FI(I)))
      ALFA=ALFAP*H(I)
      BETA=BETAP*H(I)
      QA=Q(I)/6
      TAU(I)=TAUP*TAU(I)
      BST=BETA*AST/2.
      BCT=BETA*ACT/2.
      AST=ALFA*AST
      ACT=ALFA*ACT/2.
      EIEPSZI=CEXP(CMPLX(0.,EPSZI(I)))
      S2ET=SIN(ETETA(I))**2
      CEF=COS(EFI(I))
      SEF=SIN(EFI(I))
      CET=COS(ETETA(I))
C
C     ********************* GROUND STATE MATRIX (B) *******************
C
      B(1,1,I)=CMPLX(BCT-CHEM(I),0.)
      B(1,2,I)=-BST*EIFI
      B(2,1,I)=CONJG(B(1,2,I))
      B(2,2,I)=CMPLX(-BCT-CHEM(I),0.)
C
C     ********************* EXCITED STATE MATRIX (A) ******************
C
      A(1,4,I)=(0.,0.)
      A(4,1,I)=(0.,0.)
      A(1,3,I)=QA*SQ32*(-3.*S2ET+ETA(I)*(CMPLX(CEF,SEF*CET)**2
     1                 -CMPLX(-SEF,CEF*CET)**2))*EIEPSZI**2
      A(2,4,I)=A(1,3,I)
      A(3,1,I)=CONJG(A(1,3,I))
      A(4,2,I)=A(3,1,I)
      A(2,3,I)=-AST*EIFI
      A(3,2,I)=CONJG(A(2,3,I))
      A(1,2,I)=SQ32*A(2,3,I)+QA*SQ3*(1.5*IM*SIN(2.*ETETA(I))-ETA(I)*
     1(CMPLX(CEF,SEF*CET)*SEF-CMPLX(-SEF,CEF*CET)*CEF)*SIN(ETETA(I)))*
     2EIEPSZI
      A(3,4,I)=-A(1,2,I)+SQ3*A(2,3,I)
      A(2,1,I)=CONJG(A(1,2,I))
      A(4,3,I)=CONJG(A(3,4,I))
      FF=6.*ETA(I)*S2ET*COS(2.*EFI(I))
      Z=CMPLX(0.25*QA*(27.*COS(ETETA(I))**2+9.*S2ET-15.-FF),0.)
      A(1,1,I)=Z+3.*ACT
      A(4,4,I)=Z-3.*ACT
      Z=CMPLX(0.25*QA*(3.*COS(ETETA(I))**2+21.*S2ET-15.+FF),0.)
      A(2,2,I)=Z+ACT
  20  A(3,3,I)=Z-ACT
C
C     *************** DETERMINE STOCHASTIC MATRIX ELEMENTS (W) ********
C
      W(1,2)=CMPLX(1./TAU(1),0.)
      W(2,1)=CMPLX(1./TAU(2),0.)
      W(1,1)=-W(1,2)
      W(2,2)=-W(2,1)
C
C     ****************** OBTAIN RELAXATION MATRIX (SP=SPU) ************
C
      DO 25 I=1,16
      DO 25 J=1,16
  25  SPU(I,J)=(0.,0.)
      DO 45 I=1,2
      DO 45 J=1,2
      DO 45 K=1,4
      DO 30 KV=1,4
  30  SP(I,J,K,I,J,KV)=IM*A(KV,K,I)
      DO 35 JV=1,2
  35  SP(I,J,K,I,JV,K)=SP(I,J,K,I,JV,K)-IM*B(J,JV,I)
      DO 40 IV=1,2
  40  SP(I,J,K,IV,J,K)=SP(I,J,K,IV,J,K)-W(I,IV)
  45  SP(I,J,K,I,J,K)=SP(I,J,K,I,J,K)+CMPLX(GAMMAFEL,0.)
C
C     *************************** SPECTRUM ****************************
C
      WRITE(*,'(A)') '+Eigensystem analysis...'
      CALL CMATEIG(SPU,SV,SE)
      WRITE(*,'(A)') '+Matrix inversion...'
      CALL CMATINV(SV,ISV)
C           ******************** CHANNEL LOOP *******************
      WRITE(*,46) NCHAN
  46  FORMAT(/,1X,I3,'-channels spectrum computing in progress....',/)
      DO 90 L=1,NCHAN
      WRITE(*,47) L
  47  FORMAT('+',55X,'channel = ',I3)
      OMEGA(L)=REF+(L-1)*DELTA
      DO 50 I=1,16
  50  SEO(I)=1/(SE(I)+CMPLX(0.,-OMEGA(L)))
      DO 55 I=1,16
      DO 55 J=1,16
      D(I,J)=(0.,0.)
  55  D(I,J)=SV(I,J)*SEO(J)
      DO 65 I=1,16
      DO 65 J=1,16
      SPU(I,J)=(0.,0.)
      TEMP=(0.,0.)
      DO 60 K=1,16
  60  TEMP=TEMP+D(I,K)*ISV(K,J)
  65  SPU(I,J)=TEMP
C
C     ******* AVERAGE OVER THE STATISTICAL INDICES TO OBTAIN THE ******
C     ***************** QUANTUM-MECHANICAL MATRIX (SPQM) **************
C
      DO 70 I=1,2
  70  PW(I)=CMPLX(TAU(I)/(TAU(1)+TAU(2)),0.)
      DO 75 J=1,2
      DO 75 K=1,4
      DO 75 JV=1,2
      DO 75 KV=1,4
  75  SPQM(J,K,JV,KV)=(0.,0.)
      DO 80 J=1,2
      DO 80 K=1,4
      DO 80 JV=1,2
      DO 80 KV=1,4
      SPQ=(0.,0)
      DO 77 N=1,2
      DO 77 M=1,2
  77  SPQ=SPQ+PW(N)*SP(N,J,K,M,JV,KV)
  80  SPQM(J,K,JV,KV)=SPQ
C
C     ******************** OBTAIN THE SPECTRUM (SS) *******************
C
      Z=(0.,0.)
      DO 85 J=1,2
      DO 85 K=1,4
      Z01=AP(1,J,K)
      Z02=AP(2,J,K)
      DO 85 JV=1,2
      DO 85 KV=1,4
      IF(IABS(K-J-1).EQ.2.OR.IABS(KV-JV-1).EQ.2) GO TO 85
      Z=Z+(Z01*CONJG(AP(1,JV,KV))+Z02*CONJG(AP(2,JV,KV)))*
     1SPQM(J,K,JV,KV)
  85  CONTINUE
  90  SS(L)=REAL(Z)/GAMMAFEL
      WRITE(*,'(A)') '+End of spectrum computing'
C
C     ************ CONVOLUTION ************
C
      WRITE(*,680)
  680 FORMAT(/' Convolution computing')
      RGAM22=REFGAM*REFGAM/4.
      DO 120 I=1,256
  120 SSC(I)=0.
      DO 150 L=1,NCHAN
      WRITE(*,653) L
  653 FORMAT('+Channel :',I3)
      E0=REF+(L-1)*DELTA
      DO 150 I=1,NCHAN
      EI=REF+(I-1)*DELTA
  150 SSC(L)=SSC(L)+(RGAM22/(RGAM22+(EI-E0)*(EI-E0)))*SS(I)
C     ************ AMPLITUDES DOWN *********
      SMAX=SSC(1)
      DO 160 I=2,NCHAN
      IF(SMAX.GT.SSC(I)) GOTO 160
      SMAX=SSC(I)
  160 CONTINUE
      DO 180 I=1,NCHAN
  180 SS(I)=(SMAX+0.00001-SSC(I))*100000.
C
C     ************************* DATA OUTPUT ***************************
C
      WRITE(6,642)
      WRITE(*,640) STIT,(SS(I),I=1,NCHAN)
      WRITE(6,640) STIT,(SS(I),I=1,NCHAN)
 640  FORMAT(//,1X,A72/,' RESULT SPECTRUM', //(8F10.0))
      CALL SPUTUNF(SS,NCHAN,OSPFN,STIT)
C
C     *************************    EXIT     ***************************
C
      WRITE(*,650)
      WRITE(6,642)
 642  FORMAT('1')
 650  FORMAT(//)
 100  STOP
      END
C
C     #################################################################
C     ************************ SUBROUTINE MENU ************************
C
      SUBROUTINE MENU(ANS,DATAFN,ISPFN,OSPFN,STIT)
      CHARACTER*72 STIT
      CHARACTER*12 DATAFN,ISPFN,OSPFN
      CHARACTER*1 ANS,MANS
      CALL CLS
      ANS='N'
      MANS='N'
   2  WRITE(*,5)
   5  FORMAT(////,26X,'R E L A X - M E N U',/,26X,
     1'-------------------',/////,20X,
     2'[1]  SPECTRUM TITLE SPECIFICATION',/,20X,
     3'[2]  INPUT DATA-FILE SPECIFICATION',/,20X,
     4'[3]  OUTPUT SPECTRUM-FILE SPECIFICATION',/,20X,
     5'[4]  COMPUTING SPECTRUM',/,20X,
     6'[5]  QUIT',////,2X,'enter your choise: ')
      READ(*,'(A)') MANS
      IF ((ICHAR(MANS).LT.49).OR.(ICHAR(MANS).GT.53)) GO TO 50
      IF (MANS.EQ.'1') GO TO 100
      IF (MANS.EQ.'2') GO TO 200
      IF (MANS.EQ.'3') GO TO 300
      IF (MANS.EQ.'4') GO TO 400
      IF (MANS.EQ.'5') GO TO 350
  50  WRITE (*,55)
  55  FORMAT(//,15X,' ####  NOT VALID CHOISE  ####',//,2X,
     1'press ENTER to continue')
      READ(*,'(A)') MANS
      CALL CLS
      GO TO 2
 100  WRITE(*,'(A)') ' ENTER SPECTRUM TITLE <max. 72 char> :'
      READ(*,'(A)') STIT
      GO TO 2
 200  WRITE(*,'(A)') ' ENTER INPUT DATA-FILE NAME :'
      READ(*,'(A)') DATAFN
      GO TO 2
 300  WRITE(*,'(A)') ' ENTER OUTPUT SPECTRUM FILE NAME :'
      READ(*,'(A)') OSPFN
      GO TO 2
 350  WRITE(*,'(A)') ' Quit RELAX and return to DOS ?  [Y/N] :'
      READ(*,'(A)') ANS
      IF ((ANS.EQ.'Y').OR.(ANS.EQ.'y')) GO TO 420
      GO TO 2
 400  CALL CLS
      WRITE(*,603) STIT,DATAFN,OSPFN
  603 FORMAT(//,' SPECTRUM TITLE : ',A72,/,
     1' INPUT DATA-FILE : ',A12,//,
     2' OUTPUT SPECTRUM FILE : ',A12,///,10X,
     3'LOOK AT THE INPUT OPTIONS !       O.K. ?   [Y/N]')
      READ(*,'(A)') MANS
      IF ((MANS.EQ.'Y').OR.(MANS.EQ.'y')) GO TO 420
      GO TO 2
 420  CALL CLS
      RETURN
      END
C
C     *********************** SUBROUTINE CMATEIG **********************
C
      SUBROUTINE CMATEIG(SPU,SV,SE)
      COMPLEX SPU(16,16),SV(16,16),SE(16),ZN
      REAL WK(600)
      N=16
      IA=16
      IZ=16
      IJOB=2
      CALL EIGCC (SPU,N,IA,IJOB,SE,SV,IZ,WK,IER)
      WRITE(*,5) WK(1)
      WRITE(6,5) WK(1)
   5  FORMAT(' Performance index of the eig.sys. analysis:',F10.5/)
C     NORMALIZE EIGENVECTORS
C     DO 10 J=1,N
C     ZN=SV(1,J)
C     DO 10 I=1,N
C     SV(I,J)=SV(I,J)/ZN
C 10  CONTINUE
      RETURN
      END
C
C     ********************* SUBROUTINE CMATINV ************************
C
      SUBROUTINE CMATINV(SV,ISV)
      COMPLEX SV(16,16),ISV(16,16)
      REAL A(32,32),B(32,16),WKAREA(1200)
      M=16
      N=32
      IA=32
      IDGT=0
      DO 5 I=1,16
      DO 5 J=1,16
      A(I,J)=REAL(SV(I,J))
   5  A(I+16,J+16)=A(I,J)
      DO 10 I=1,16
      DO 10 J=1,16
      A(I+16,J)=AIMAG(SV(I,J))
  10  A(I,J+16)=-A(I+16,J)
      DO 15 I=1,32
      DO 15 J=1,16
  15  B(I,J)=0.0
      DO 20 I=1,16
  20  B(I,I)=1.0
      CALL LEQT2F (A,M,N,IA,B,IDGT,WKAREA,IER)
      DO 25 I=1,16
      DO 25 J=1,16
  25  ISV(I,J)=CMPLX(B(I,J),B(I+16,J))
      RETURN
      END
C
C     ******************** SUBROUTINE SPUTUNF *************************
C
      SUBROUTINE SPUTUNF(SS,NCHAN,OSPFN,STIT)
      DIMENSION SS(256)
      CHARACTER*12 OSPFN
      CHARACTER*72 STIT
      OPEN(3,FILE='RELAX2X.UNF',STATUS='NEW',FORM='UNFORMATTED')
      WRITE(3) NCHAN,STIT
      WRITE(3) (SS(I),I=1,NCHAN)
      CLOSE(3,STATUS='KEEP')
      WRITE(*,641)
      WRITE(6,641)
 641  FORMAT(' Unformatted file of RELAX2X.UNF was created'/)
      OPEN(9,FILE=OSPFN,STATUS='NEW',FORM='FORMATTED')
      WRITE(9,20) STIT
  20  FORMAT(A72)
      WRITE(9,30) (SS(I),I=1,NCHAN)
  30  FORMAT(F12.0)
      CLOSE(9,STATUS='KEEP')
      WRITE(*,642) OSPFN
      WRITE(6,642) OSPFN
  642 FORMAT(' The result spectrum was written to file of: ',A12/)
      RETURN
      END
C
C     ************************ SUBROUTINE CLS *************************
C
      SUBROUTINE CLS
      CHARACTER*1 ESC
      ESC=CHAR(27)
      WRITE(*,1)ESC
   1  FORMAT(1X,A1,'[2J')
      RETURN
      END
