c **Simulates polycrystalline e.s.r spectra - includes
c ** g and a tensors with coincident principal axes - 
c **ligand hyperfine splittings due to one type of 
c **nuclei can be included in first order - the ligand
c **tensors are assumed to be axial with principal
c **direction specified in the input. Second order 
c **perturbation equations are used to calculate line
c **position. 3-point quadrature formula is used for
c **integration. Lorentzian or Gaussian shape 
c **with anisotropic M(I)-dependent width is allowed. 
c ** M(I)-dependent width is used. 
c **M. V. Rajasekharan, Madras 1977,
c **modified, Hyderabad, 1987, 1992, 1994.
c **See PSPM.DOC for input info.
c
c ********FILES
c                 UNIT 5    input data
c                 *         output data
c                 UNIT 8    x,y data for plotting using
c                           ESRS.BAS (screen) or
c                           ESRP.BAS (plotter) 
C------------------------------------------------------
C **(1) arrays and constants:
C
      dimension xx(3),a(3)
      dimension se(16),sn(16),f3(16),f4(16),f5(16),
     1f6(16),f3r(16),f4r(16),f5r(16),f6r(16),tmp(16),
     1tmm(16),tmh(16),tml(16),lg(120),le(120)
      dimension dlx(10),dly(10),dlz(10),alig(10),
     1blig(10),alg(10),nlsub(10)
      CHARACTER*10 TITLE,FPLOT
      CHARACTER*2 SB,INS
      DIMENSION SB(26)
      COMMON X(1200),Y(1200)
      data ngl,xx,a/3,0.7745966692,0.0,-0.7745966692,
     10.5555555556,0.8888888889,0.5555555556/
	DATA SB/'SC','GX','GY','GZ','AX','AY','AZ','QP',
     1'QD','GN','FR','TS','PS','SF','EF','WI','SL',
     1'LP','TE','PZ','WX','WY', 'WZ','AF','GM','LS'/
      data beta,betan/4.668593e-05,2.5426742e-08/
      conv=3.1415927/180.0
      cong=714.5045*29.980
      PFAC=10.0
      PRN=1020.0
      N=1021
C
C **(2) set default values:
C 
      NEBP=-2
      TRA=0
      PRA=0
      NPCHK=1000
	MULTC=1
	GX=2.
	GY=2.
	GZ=2.
	AX=0.1E-7
	AY=0.1E-7
	AZ=0.1E-7
	QP=0.
	QDP=0.
	GN=0.
	FR=0.3
	TR=90.
	PR=180.
        TE=90.
        PE=180.
	HST=2500.
	HSTC=2500.
	HED=3500.
	HEDC=3500.
	WX=5.
        WY=5.
        WZ=5.
        ALFA=0
        GAMA=0
	MULTL=1
	NLIG=0
	L=1
C
C **(3) read input from unit 5:
C
	READ(5,300)TITLE
        READ(5,300)FPLOT
 300    FORMAT(A10)
 315    READ(5,311,ERR=318)INS,V
 311    FORMAT(A2,F12.6)
	IF(INS.EQ.'GO')GO TO 309
        IF(INS.EQ.'CO')GO TO 319
	IF(INS.EQ.'**')GO TO 315
	DO 312 I=1,26
	IF(INS.NE.SB(I))GO TO 312
	J=I
	GO TO 314
 312    CONTINUE
 318    WRITE(*,313)INS
 313    FORMAT(5X,'BAD INSTRUCTION',A2)
	STOP
 314    GO TO(321,322,323,324,325,326,327,328,329,330,
     1  331,332,333,334,335,336,337,338,339,341,
     2  342,343,344,345,346,347),J
 321    MULTC=(2.*V+1)
	GO TO 315
 322    GX=V
	GO TO 315
 323    GY=V
	GO TO 315
 324    GZ=V
	GO TO 315
 325    AX=V/10000.
	GO TO 315
 326    AY=V/10000.
	GO TO 315
 327    AZ=V/10000.
	GO TO 315
 328    QP=V/10000.
	GO TO 315
 329    QDP=V/10000.
	GO TO 315
 330    GN=V
	GO TO 315
 331    FR=V
	GO TO 315
 332    TR=V
	GO TO 315
 333    PR=V
	GO TO 315
 334    HST=V
	HSTC=V
	GO TO 315
 335    HED=V
	HEDC=V
	GO TO 315
 336    WID=V
        WX=V
        WY=V
        WZ=V 
	GO TO 315
 337    MULTL=2.*V+1
	GO TO 315
 338    READ(5,316,ERR=318)DLX(L),DLY(L),DLZ(L),
     1  ALIG(L),BLIG(L)
        ALIG(L)=ALIG(L)/10000.0
        BLIG(L)=BLIG(L)/10000.0
        NLIG=NLIG+1
        L=L+1
        GO TO 315
 339    TE=V
        GO TO 315
 341    NPCHK=V
        GO TO 315
 342    WX=V
        GO TO 315
 343    WY=V
        GO TO 315
 344    WZ=V
        GO TO 315
 345    ALFA=V
        GO TO 315
 346    GAMA=V/10000.
        GO TO 315         
 347    LS=V
        GO TO 315	     
 316    FORMAT(F12.6)
 319    NEBP=0
 309    CONTINUE
        ALFP=ALFA/BETA
        GAMP=GAMA/BETA	   
        ALFN=ALFP*FR
      gbn = gn*betan
      px = qdp - qp/3.
      py = -qdp - qp/3.
      pz = 2.*qp/3.
C
C **(4) write input to '*'(screen):
C
      write(*,300) title
   52 format(2X,'2I+1 central atom :',I2,2X,
     1', ligand :',I2,2x,'No. of ligands =',I2)
      write(*,52) multc,multl,nlig
      write(*,40) gx,gy,gz,ax,ay,az,qp,qdp,gn,fr
   40 format(/1x,'gx = ',f10.5,5x,'gy = ',f10.5,5x,
     1 'gz = ',f10.5/1x,'ax = ',f10.5,5x,'ay = ',f10.5,
     2 5x,'az = ',f10.5/1x,'qp = ',f10.5,5x,'qdp = ',
     3 f10.5,5x,'gn = ',f10.5/1x,'frequency = ',f10.5,
     4 ' cm inverse')
      if(ax.eq.0.) ax=0.1e-7
      if(ay.eq.0.) ay=0.1e-7
      if(az.eq.0.) az=0.1e-7
      if(nlig.le.0) go to 69
      do 61 i=1,nlig
      write(*,62) i,multl,dlx(i),dly(i),dlz(i),
     1alig(i),blig(i)
      alig(i)=alig(i)*alig(i)
      blig(i)=blig(i)*blig(i)
   61 continue
   62 format(1X,' l',i2,'  mult',i2,' /cosines',3f8.5,
     1 '/  a',f8.5,' b',f8.5)
   69 continue
      write(*,30)tr,pr
   30 FORMAT(1X,'THETA STEP = ',F10.6,2X,
     1 'PHI STEP = ',F10.6)
C
C **(5)branch to read previously accumulated data if 
C **'CO' instruction was given(ie., NEBP=0), 
C 
   53 IF(NEBP) 206,207,207
C **else initialise x,y data:
C  
  206 OPEN(8,FILE=FPLOT,STATUS='NEW')
      gap = (hedc-hstc)/prn
      hy = hstc
      do 208 i=1,n
      y(i)=hy
      hy = hy + gap
  208 x(i) = 0.0
      go to 209
  207 CONTINUE
      CALL PLTF(FPLOT,0,N,HST,HED,TRA)
      GAP=(HED-HST)/PRN
  209 continue
      KT=(TE-TRA)/TR
      kp = 180.0/pr 
      NPCHK=(NPCHK/KP)*KP
      IF (NPCHK.EQ.0) NPCHK=KP	
      write(*,45) hst,hed,WX,WY,WZ,ALFA,GAMA,KT,KP,NPCHK
   45 format(2x,'starting field = ',F10.3,2X,
     1 'ending field = ',F10.3/' WX = ',F10.3, ' WY = ',
     2 F10.3,' WZ = ',F10.3/' ALFA = ',E11.4,'  GAMA = ',
     3 E11.4/' KT = ',I3,'  KP = ',I3,' NPCHK = ',I4//)
      GO TO(46,47),LS
   46 WRITE(*,48)
      WFAC=1.
      GOTO 59
   47 WRITE(*,49)
      WFAC=1.2011
      PFAC=PFAC/3.3333
   48 FORMAT('  Lorentzian shape'//)
   49 FORMAT('  Gaussian shape'//)
C  
C **(6) generate spin functions: they are ordered as
C ** |1/2 I>,|1/2 I-1>,..|1/2 -I>, |-1/2 -I>,
C ** |-1/2, -I+1>,..|-1/2,I>:
C   
   59 mc = multc
      spnc = multc - 1
      spnc = spnc/2.0
      svcc = spnc*(spnc + 1.0)
      spnl = multl - 1
      spnl = spnl/2.0
      svcl = spnl*(spnl + 1.0)
      nlprod = 1
      IF(NLIG.GT.0) THEN
      DO 58 I=1,NLIG
   58 NLPROD=NLPROD*MULTL
      ENDIF
      do 51 i = 1,mc
      zz = i - 1
      sni = spnc - zz
      se(i) = 0.5
      se(i + mc) = -0.5
      sn(i) = sni
   51 sn(i + mc) = -sni
C
C **(7) generate I-, I+ multipliers etc.:
C
      mc2 = mc*2
      do 20 i = 1,mc2
      sni = sn(i)
      sei = se(i)
      f3(i) = svcc - sni*(sni + 1.0)
      f4(i) = svcc - sni*(sni - 1.0)
      f5(i) = f3(i)*(svcc - (sni + 1.0)*(sni + 2.0))
      f6(i) = f4(i)*(svcc - (sni - 1.0)*(sni - 2.0))
      tmp(i) = 3.0*(1.0 + sni)
      tmm(i) = 3.0*(1.0 - sni)
      tmh(i) = 3.0*(sni + 0.5)
      tml(i) = 3.0*(sni - 0.5)
   20 continue
      do 25 i = 1,mc
      lg(i) = mc - i + 1
   25 le(i) = mc + i
C
C **(8) start spectral calculation:
C
      npnt= 0
      nchk = npchk
      tra = tra - tr
      ppra = pra - pr
C
C **(9) theta sum:
C 
      do 100 it = 1,kt
      tra = tra + tr
      pra = ppra
C
C **(10) phi sum:
C
      WRITE(*,454)NPNT
  454 FORMAT(1H+I4)
      if(npnt.lt.nchk) go to 451
      nchk= nchk+npchk
      write(*,452) npnt,it,kt,jt,kp
  452 format (1x,i6,'points done ','ntha = ',i3,
     1 ' (',i3,')',' nphi = ',i3,' (',i3,')')
  453 format (i2)
      CALL PLTF(FPLOT,1,N,HST,HED,TRA)
  451 continue
      do 101 jt = 1,kp
      pra = pra + pr
C
C **(11) three point quadrature (NGL=3)
C 
C **     over theta:
C
      do 50 i = 1,ngl
      thetd = tra + (xx(i) + 1.0)*tr/2.0
      theta = conv*thetd
      sta = sin(theta)
      cta = cos(theta)
C 
C **     over phi:
C 
      do 50 j = 1,ngl
      pfid = pra + (xx(j) + 1.0)*pr/2.0
      pfi = conv*pfid
      spf = sin(pfi)
      cpf = cos(pfi)
C
C **(12) calculate g:
C	
      hx = sta*cpf
      hy = sta*spf
      hz = cta
      g = sqrt(hx*hx*gx*gx + hy*hy*gy*gy + hz*hz*gz*gz)
      gb = g*beta
C
C **(13) calculate A(ligand) for all ligands:
C
      if(nlig.le.0) go to 63
      do 64 ilg = 1,nlig
      costl = (dlx(ilg)*hx + dly(ilg)*hy + dlz(ilg)*hz)
      costl = costl*costl
      sintl = 1.0 - costl
      atlg = sqrt(costl*alig(ilg) + sintl*blig(ilg))
      alg(ilg) = atlg/gb
   64 continue
C
C **(14) calculate components of spin rotations,
C **l1, l2, l3, m1,m2, n1, n2, n3, L1, L2, L3, 
C **M1, M2, N1, N2, N3 as per notes:
C
   63 continue
      if(hz - 1.0) 13,11,13
   13 a11 = hx*gx/g
      a12 = hy*gy/g
      a13 = hz*gz/g
      dena = sqrt(a11*a11 + a12*a12)
      a21 = -a12/dena
      a22 = a11/dena
      a31 = a11*a13/dena
      a32 = a12*a13/dena
      a33 = -dena
C
C ** first order metal h.f, and anisotropic width:
C
      AA=SQRT((A11*AX)**2+(A12*AY)**2+(A13*AZ)**2)
      WIDT=((A11*WX*GX)**2+(A12*WY*GY)**2+
     1 (A13*WZ*GZ)**2)/(G*G)  
C
      ga=g*aa
      b11 = hx*gx*ax/ga
      b12 = hy*gy*ay/ga
      b13 = hz*gz*az/ga
      denb = sqrt(b11*b11 + b12*b12)
      b21 = -b12/denb
      b22 = b11/denb
      b31 = b11*b13/denb
      b32 = b12*b13/denb
      b33 = -denb
C
C ** anisotropic transition probability
C
      tsa = 1.0 - hz*hz
      tca = hz*hz
      psa = hy*hy/tsa
      pca = hx*hx/tsa
      pmi = (gx*gx*gy*gy*tsa + gy*gy*gz*gz*
     1 (psa + tca*pca) + gx*gx*gz*gz*
     2 (pca + tca*psa))/(g*g)
C
      go to 12
C	
C **special case, H || Z  
C 
   11 a13 = 1.
      a22 = 1.
      a31 = 1.
      b13 = 1.
      b22 = 1.
      b31 = 1.
      a11 = 0.
      a21 = 0.
      a32 = 0.
      a12 = 0.
      a33 = 0.
      b11 = 0.
      b21 = 0.
      b32 = 0.
      b12 = 0.
      b33 = 0.
      pmi = (gx*gx + gy*gy)/2.0
      AA=AZ
      WIDT=WZ*WZ
C
C **weighted intensity factor and factors for M(I)
C **dependent widths:	
C
   12 continue
      pin = pmi*a(i)*a(j)*sta
      HFWA=ALFN/(G*G)
      HFWB=GAMP/G - AA*ALFP/(G*G)
C
C **various other orientation dependent quantities:
C
      t33 = aa
      q33 = b11*b11*px + b12*b12*py + b13*b13*pz
      z3 = hx*b11 + hy*b12 + hz*b13
      t11 = a31*b31*ax + a32*b32*ay + a33*b33*az
      t22 = a21*b21*ax + a22*b22*ay
      t12 = a31*b21*ax + a32*b22*ay
      t21 = a21*b31*ax + a22*b32*ay
      t13 = a31*b11*ax + a32*b12*ay + a33*b13*az
      t23 = a21*b11*ax + a22*b12*ay
      q11 = b31*b31*px + b32*b32*py + b33*b33*pz
      q22 = b21*b21*px + b22*b22*py
      q12 = b31*b21*px + b32*b22*py
      q13 = b31*b11*px + b32*b12*py + b33*b13*pz
      q23 = b21*b11*px + b22*b12*py
      z1 = hx*b31 + hy*b32 + hz*b33
      z2 = hx*b21 + hy*b22
      t33s = 0.5*t33
      qms = ((q11 - q22)**2 + 4.0*q12*q12)/32.0
      tps = t13*t13 + t23*t23
      tp = (t11 + t22)**2
      tpp = (t12 + t21)**2
      tm = (t11 - t22)**2
      tpm = (t21 - t12)**2
C
C **sum over metal h.f components, delta-M=0:
C
      do 60 ki = 1,mc
      k = le(ki)
      snk = sn(k)
C
C **M(I) dependent width
C 
      WIDH=HFWA+HFWB*SNK
      WID=WFAC*SQRT(WIDT+WIDH*WIDH)
      WID2=WID*WID
      PINL=PIN*WID
      PING=PIN/(WID*WID2)
      WNW=WID*PFAC
C
C **first order field	
C
      h1s = (fr - t33*snk)/gb
C
C **second order corretion
C
      c2i = 2.0*qms*(f5(k) - f6(k))/t33s
      gbh = gb*h1s
      gbnh = gb*h1s*z3
      gbns = gbn*h1s*z1
      gbnt = gbn*h1s*z2
      t1 = f4(k)*((2.*q13*snk - q13 - gbns)**2 +
     1 (2.*q23*snk - q23 - gbnt)**2)/(2.*t33s)
      t2 = f3(k)*((2.*q13*snk + q13 - gbns)**2 +
     1 (2.*q23*snk + q23 - gbnt)**2)/(2.*t33s)
      t3 = snk*snk*tps/(2.*(gbh + t33*snk))
      t4 = (tm + tp + tpm + tpp)*(f3(k)/
     1 (16.*(gbh + t33*(snk + 0.5)))
     2 + f4(k)/(16.*(gbh + t33*(snk - 0.5))))
      dh2 = (c2i - t1 + t2 - t3 - t4)/gb
      h2c = h1s + dh2
C
C **add ligand h.f splittings:
C    
      do 65 jlg = 1,nlprod
      h2 = h2c
      if(nlig.le.0) go to 68
      mlg = jlg
      jpr = nlprod
      do 66 klg=1,nlig
      jpr = jpr/multl
      nlsub(klg) = (mlg-1)/jpr + 1
      mlg = mod(mlg-1,jpr) + 1
   66 continue
      do 67 klg = 1,nlig
C
C **final field value of a line in the spectrum
C
      h2 = h2 + alg(klg)*(nlsub(klg) - spnl - 1.)
   67 continue
   68 continue
C 
C **field range for computing the signal with width   
C 
      npl = 1. + (h2 - wnw - hstc)/gap
      npr = 1. + (h2 + wnw - hstc)/gap
      IF(NPL.LT.1) THEN
	NPL = 1
	GOTO 56
	ENDIF
	IF(NPR.GT.N) THEN
	NPR = N
	GOTO 56
	ENDIF
C	
C **computation of a line within field limits
C
	NPM = 0.51 + FLOAT(NPR-NPL)/2.
	DO 54 IPT=1,NPM
	DHF = Y(NPL) - H2
      IF(LS.EQ.2) GOTO 71
C **Lorentzian shape    
      TRMZ = PINL*DHF/(WID2 + DHF*DHF)**2
      GOTO 55
C **Gaussian shape    
   71 TRMZ = PING*DHF*EXP(-DHF*DHF/WID2)
   55 X(NPL)= X(NPL) + TRMZ
      X(NPR)= X(NPR) - TRMZ
	NPL = NPL + 1
	NPR = NPR - 1
   54 CONTINUE      
      GOTO 65
C
C **computation of a line cutting at a border    
C
   56 do 70 ipt = npl,npr
      dhf = y(ipt) - h2
      IF (LS.EQ.2) GOTO 72
C **Lorentzian shape    
	TRMZ = PINL*DHF/(WID2 + DHF*DHF)**2
      GOTO 70
C **Gaussian shape    
   72 TRMZ = PING*DHF*EXP(-DHF*DHF/WID2)
   70 x(ipt) = x(ipt) + TRMZ
C
C **end of a line   
   65 continue
C **ligand splittings end   
   60 continue
C **central atom splittings end   
   50 continue
C **quadrature in theta and phi end   
  101 continue
C **phi sum end  
      NPNT=NPNT+KP
  100 continue
C **theta sum end  
  102 continue
C 
C **write x,y data to file:
C
      TRA=TRA+TR
      CALL PLTF(FPLOT,1,N,HST,HED,TRA)
      STOP
      END
C
      SUBROUTINE PLTF(F,L,N,S,E,T)
      CHARACTER*10 F
      COMMON X(1200),Y(1200)
      OPEN(8,FILE=F)
      REWIND 8
C
C **if L=0 reads previous data and returns:
C	
      IF(L.EQ.0) GO TO 10
C
C **else writes present data and returns:
C **find highest and lowest X values(X contains
C **intensities and Y the field values)
C
      XL=X(1)
      XH=XL
      DO 409 I=1,N
      IF(XL-X(I))410,410,411
  411 XL=X(I)
  410 IF(XH-X(I))412,409,409
  412 XH=X(I)
  409 CONTINUE
C
C **write the no. of points; starting and ending 
C **fields, X-low and X-high; current theta and phi;
C **X,Y data:
C     DUMY=0
      WRITE(8,455) N
      WRITE(8,456) S,E,XL,XH
      WRITE(8,466) T,DUMY
      WRITE(8,407) (Y(I),X(I),I=1,N)
      RETURN
C
C **read no. of points; starting and ending fields;
C **theta and phi; X, Y data: (NOTE: S and E read  
C **here will override the input by 'SF' and 'EF') 
C	
   10 READ(8,455) N
      READ(8,456) S,E,DUMC,DUMD
      READ(8,466) T,DUMY
      READ(8,407) (Y(I),X(I),I=1,N)
      RETURN
  455 FORMAT(I5)
  456 FORMAT(4(E11.4,1X))
  466 FORMAT(2(E11.4,1X))
  407 FORMAT(3(2X,E11.4,1X,E11.4))
      END