text:
c######################################################################
c----------------------------------------------------------------------
      subroutine easypbe(up,agrup,delgrup,uplap,dn,agrdn,delgrdn,dnlap,
     1           agr,delgr,lcor,lpot,
     1           exlsd,vxuplsd,vxdnlsd,eclsd,vcuplsd,vcdnlsd,
     1           expw91,vxuppw91,vxdnpw91,ecpw91,vcuppw91,vcdnpw91,
     1           expbe,vxuppbe,vxdnpbe,ecpbe,vcuppbe,vcdnpbe)
c----------------------------------------------------------------------
c----------------------------------------------------------------------
c EASYPBE is a driver for the PBE subroutines, using simple inputs
c K. Burke, May 14, 1996.
c inputs: up=up density
c	: agrup=|grad up|
c	: delgrup=(grad up).(grad |grad up|)
c	: uplap=grad^2 up=Laplacian of up
c	: dn,agrdn,delgrdn,dnlap=corresponding down quantities
c	: agr=|grad rho|
c	: delgr=(grad rho).(grad |grad rho|)
c	: lcor=flag to do correlation(=0=>don't)
c	: lpot=flag to do potential(=0=>don't)
c outputs: exlsd=LSD exchange energy density, so that
c		ExLSD=int d^3r rho(r) exlsd(r)
c	 : vxuplsd=up LSD exchange potential
c	 : vxdnlsd=down LSD exchange potential
c        : exclsd=LSD exchange-correlation energy density
c	 : vxcuplsd=up LSD exchange-correlation potential
c	 : vxcdnlsd=down LSD exchange-correlation potential
c        : expw91,vxuppw91,vxdnpw91,ecpw91,etc.=PW91 quantities
c        : expbe,vxuppbe,vxdnpbe,ecpbe,etc.=PBE quantities
c----------------------------------------------------------------------
c----------------------------------------------------------------------
c needed constants:
c pi32=3 pi**2
c alpha=(9pi/4)**thrd
      implicit real*8(a-h,o-z)
      parameter(thrd=1.d0/3.d0,thrd2=2.d0*thrd)
      parameter(pi32=29.608813203268075856503472999628d0)
      parameter(pi=3.1415926535897932384626433832795d0)
      parameter(alpha=1.91915829267751300662482032624669d0)
c----------------------------------------------------------------------
c----------------------------------------------------------------------
c PBE exchange
c use  Ex[up,dn]=0.5*(Ex[2*up]+Ex[2*dn]) (i.e., exact spin-scaling)
c do up exchange
c fk=local Fermi wavevector for 2*up=(3 pi^2 (2up))^(1/3)
c s=dimensionless density gradient=|grad rho|/ (2*fk*rho)_(rho=2*up)
c u=delgrad/(rho^2*(2*fk)**3)_(rho=2*up)
c v=Laplacian/(rho*(2*fk)**2)_(rho=2*up)
      rho2=2.d0*up
      if(rho2.gt.1d-18)then
        fk=(pi32*rho2)**thrd
        s=2.d0*agrup/(2.d0*fk*rho2)
        u=4.d0*delgrup/(rho2*rho2*(2.d0*fk)**3)
        v=2.d0*uplap/(rho2*(2.d0*fk)**2)
        call exchpbe(rho2,s,u,v,0,lpot,exuplsd,vxuplsd)
        call exchpw91(rho2,s,u,v,exuppw91,vxuppw91)
        call exchpbe(rho2,s,u,v,1,lpot,exuppbe,vxuppbe)
      else
	exuplsd=0.d0
	vxuplsd=0.d0
	exuppw91=0.d0
	vxuppw91=0.d0
	exuppbe=0.d0
	vxuppbe=0.d0
      endif
c repeat for down
      rho2=2.d0*dn
      if(rho2.gt.1d-18)then
        fk=(pi32*rho2)**thrd
        s=2.d0*agrdn/(2.d0*fk*rho2)
        u=4.d0*delgrdn/(rho2*rho2*(2.d0*fk)**3)
        v=2.d0*dnlap/(rho2*(2.d0*fk)**2)
        call exchpbe(rho2,s,u,v,0,lpot,exdnlsd,vxdnlsd)
        call exchpw91(rho2,s,u,v,exdnpw91,vxdnpw91)
        call exchpbe(rho2,s,u,v,1,lpot,exdnpbe,vxdnpbe)
      else
	exdnlsd=0.d0
	vxdnlsd=0.d0
	exdnpw91=0.d0
	vxdnpw91=0.d0
	exdnpbe=0.d0
	vxdnpbe=0.d0
      endif
10    continue
c construct total density and contribution to ex
      rho=up+dn
      exlsd=(exuplsd*up+exdnlsd*dn)/rho
      expw91=(exuppw91*up+exdnpw91*dn)/rho
      expbe=(exuppbe*up+exdnpbe*dn)/rho
      if(lcor.eq.0)return
c----------------------------------------------------------------------
c----------------------------------------------------------------------
c Now do correlation
c zet=(up-dn)/rho
c g=phi(zeta)
c rs=(3/(4pi*rho))^(1/3)=local Seitz radius=alpha/fk
c sk=Ks=Thomas-Fermi screening wavevector=sqrt(4fk/pi)
c twoksg=2*Ks*phi
c t=correlation dimensionless gradient=|grad rho|/(2*Ks*phi*rho)
c uu=delgrad/(rho^2*twoksg^3)
c rholap=Laplacian
c vv=Laplacian/(rho*twoksg^2)
c ww=(|grad up|^2-|grad dn|^2-zet*|grad rho|^2)/(rho*twoksg)^2
c ec=lsd correlation energy
c vcup=lsd up correlation potential
c vcdn=lsd down correlation potential
c h=gradient correction to correlation energy
c dvcup=gradient correction to up correlation potential
c dvcdn=gradient correction to down correlation potential
      if(rho.lt.1.d-18) then      ! modified, LWW
      vcuppbe=0.d0
      vcdnpbe=0.d0
      ecpbe=0.d0
      return
      endif
      zet=(up-dn)/rho
      g=((1.d0+zet)**thrd2+(1.d0-zet)**thrd2)/2.d0
      fk=(pi32*rho)**thrd
      rs=alpha/fk
      sk=sqrt(4.d0*fk/pi)
      twoksg=2.d0*sk*g
      t=agr/(twoksg*rho)
      uu=delgr/(rho*rho*twoksg**3)
      rholap=uplap+dnlap
      vv=rholap/(rho*twoksg**2)
      ww=(agrup**2-agrdn**2-zet*agr**2)/(rho*rho*twoksg**2)
      call CORPBE(RS,ZET,T,UU,VV,WW,1,lpot,ec,vcup,vcdn,
     1                  H,DVCUP,DVCDN)
      eclsd=ec
      ecpbe=ec+h
      vcuplsd=vcup
      vcdnlsd=vcdn
      vcuppbe=vcup+dvcup
      vcdnpbe=vcdn+dvcdn
      call CORLSD(RS,ZET,EC,VCUP,VCDN,ECRS,ECZET,ALFC)
      call CORPW91(RS,ZET,G,EC,ECRS,ECZET,T,UU,VV,WW,H,DVCUP,DVCDN)
      ecpw91=ec+h
      vcuppw91=vcup+dvcup
      vcdnpw91=vcdn+dvcdn
      return
      end
c----------------------------------------------------------------------
c######################################################################
c----------------------------------------------------------------------
      SUBROUTINE EXCHPBE(rho,S,U,V,lgga,lpot,EX,VX)
c----------------------------------------------------------------------
C  PBE EXCHANGE FOR A SPIN-UNPOLARIZED ELECTRONIC SYSTEM
c  K Burke's modification of PW91 codes, May 14, 1996
c  Modified again by K. Burke, June 29, 1996, with simpler Fx(s)
c----------------------------------------------------------------------
c----------------------------------------------------------------------
C  INPUT rho : DENSITY
C  INPUT S:  ABS(GRAD rho)/(2*KF*rho), where kf=(3 pi^2 rho)^(1/3)
C  INPUT U:  (GRAD rho)*GRAD(ABS(GRAD rho))/(rho**2 * (2*KF)**3)
C  INPUT V: (LAPLACIAN rho)/(rho*(2*KF)**2)
c   (for U,V, see PW86(24))
c  input lgga:  (=0=>don't put in gradient corrections, just LDA)
c  input lpot:  (=0=>don't get potential and don't need U and V)
C  OUTPUT:  EXCHANGE ENERGY PER ELECTRON (EX) AND POTENTIAL (VX)
c----------------------------------------------------------------------
c----------------------------------------------------------------------
c References:
c [a]J.P.~Perdew, K.~Burke, and M.~Ernzerhof, submiited to PRL, May96
c [b]J.P. Perdew and Y. Wang, Phys. Rev.  B {\bf 33},  8800  (1986);
c     {\bf 40},  3399  (1989) (E).
c----------------------------------------------------------------------
c----------------------------------------------------------------------
c Formulas:
c   	e_x[unif]=ax*rho^(4/3)  [LDA]
c ax = -0.75*(3/pi)^(1/3)
c	e_x[PBE]=e_x[unif]*FxPBE(s)
c	FxPBE(s)=1+uk-uk/(1+ul*s*s)                 [a](13)
c uk, ul defined after [a](13)
c----------------------------------------------------------------------
c----------------------------------------------------------------------
      IMPLICIT REAL*8 (A-H,O-Z)
      parameter(thrd=1.d0/3.d0,thrd4=4.d0/3.d0)
      parameter(pi=3.14159265358979323846264338327950d0)
      parameter(ax=-0.738558766382022405884230032680836d0)
      parameter(um=0.2195149727645171d0,uk=0.8040d0,ul=um/uk)
c----------------------------------------------------------------------
c----------------------------------------------------------------------
c construct LDA exchange energy density
      exunif = AX*rho**THRD
      if(lgga.eq.0)then
	ex=exunif
        vx=ex*thrd4
	return
      endif
c----------------------------------------------------------------------
c----------------------------------------------------------------------
c construct PBE enhancement factor
      S2 = S*S
      P0=1.d0+ul*S2
      FxPBE = 1d0+uk-uk/P0
      EX = exunif*FxPBE
      if(lpot.eq.0)return
c----------------------------------------------------------------------
c----------------------------------------------------------------------
C  ENERGY DONE. NOW THE POTENTIAL:
c  find first and second derivatives of Fx w.r.t s.
c  Fs=(1/s)*d FxPBE/ ds
c  Fss=d Fs/ds
      Fs=2.d0*uk*ul/(P0*P0)
      Fss=-4.d0*ul*S*Fs/P0
c----------------------------------------------------------------------
c----------------------------------------------------------------------
c calculate potential from [b](24)
      VX = exunif*(THRD4*FxPBE-(U-THRD4*S2*s)*FSS-V*FS)
      RETURN
      END
c----------------------------------------------------------------------
c######################################################################
c----------------------------------------------------------------------
      SUBROUTINE CORPBE(RS,ZET,T,UU,VV,WW,lgga,lpot,ec,vcup,vcdn,
     1                  H,DVCUP,DVCDN)
c----------------------------------------------------------------------
c  Official PBE correlation code. K. Burke, May 14, 1996.
C  INPUT: RS=SEITZ RADIUS=(3/4pi rho)^(1/3)
C       : ZET=RELATIVE SPIN POLARIZATION = (rhoup-rhodn)/rho
C       : t=ABS(GRAD rho)/(rho*2.*KS*G)  -- only needed for PBE
C       : UU=(GRAD rho)*GRAD(ABS(GRAD rho))/(rho**2 * (2*KS*G)**3)
C       : VV=(LAPLACIAN rho)/(rho * (2*KS*G)**2)
C       : WW=(GRAD rho)*(GRAD ZET)/(rho * (2*KS*G)**2
c       :  UU,VV,WW, only needed for PBE potential
c       : lgga=flag to do gga (0=>LSD only)
c       : lpot=flag to do potential (0=>energy only)
c  output: ec=lsd correlation energy from [a]
c        : vcup=lsd up correlation potential
c        : vcdn=lsd dn correlation potential
c        : h=NONLOCAL PART OF CORRELATION ENERGY PER ELECTRON
c        : dvcup=nonlocal correction to vcup
c        : dvcdn=nonlocal correction to vcdn
c----------------------------------------------------------------------
c----------------------------------------------------------------------
c References:
c [a] J.P.~Perdew, K.~Burke, and M.~Ernzerhof,
c     {\sl Generalized gradient approximation made simple}, sub.
c     to Phys. Rev.Lett. May 1996.
c [b] J. P. Perdew, K. Burke, and Y. Wang, {\sl Real-space cutoff
c     construction of a generalized gradient approximation:  The PW91
c     density functional}, submitted to Phys. Rev. B, Feb. 1996.
c [c] J. P. Perdew and Y. Wang, Phys. Rev. B {\bf 45}, 13244 (1992).
c----------------------------------------------------------------------
c----------------------------------------------------------------------
      IMPLICIT REAL*8 (A-H,O-Z)
c thrd*=various multiples of 1/3
c numbers for use in LSD energy spin-interpolation formula, [c](9).
c      GAM= 2^(4/3)-2
c      FZZ=f''(0)= 8/(9*GAM)
c numbers for construction of PBE
c      gamma=(1-log(2))/pi^2
c      bet=coefficient in gradient expansion for correlation, [a](4).
c      eta=small number to stop d phi/ dzeta from blowing up at
c          |zeta|=1.
      parameter(thrd=1.d0/3.d0,thrdm=-thrd,thrd2=2.d0*thrd)
      parameter(sixthm=thrdm/2.d0)
      parameter(thrd4=4.d0*thrd)
      parameter(GAM=0.5198420997897463295344212145565d0)
      parameter(fzz=8.d0/(9.d0*GAM))
      parameter(gamma=0.03109069086965489503494086371273d0)
      parameter(bet=0.06672455060314922d0,delt=bet/gamma)
      parameter(eta=1.d-12)
c----------------------------------------------------------------------
c----------------------------------------------------------------------
c find LSD energy contributions, using [c](10) and Table I[c].
c EU=unpolarized LSD correlation energy
c EURS=dEU/drs
c EP=fully polarized LSD correlation energy
c EPRS=dEP/drs
c ALFM=-spin stiffness, [c](3).
c ALFRSM=-dalpha/drs
c F=spin-scaling factor from [c](9).
c construct ec, using [c](8)
      rtrs=dsqrt(rs)
      CALL gcor2(0.0310907D0,0.21370D0,7.5957D0,3.5876D0,1.6382D0,
     1    0.49294D0,rtrs,EU,EURS)
      CALL gcor2(0.01554535D0,0.20548D0,14.1189D0,6.1977D0,3.3662D0,
     1    0.62517D0,rtRS,EP,EPRS)
      CALL gcor2(0.0168869D0,0.11125D0,10.357D0,3.6231D0,0.88026D0,
     1    0.49671D0,rtRS,ALFM,ALFRSM)
      ALFC = -ALFM
      Z4 = ZET**4
      F=((1.D0+ZET)**THRD4+(1.D0-ZET)**THRD4-2.D0)/GAM
      EC = EU*(1.D0-F*Z4)+EP*F*Z4-ALFM*F*(1.D0-Z4)/FZZ
c----------------------------------------------------------------------
c----------------------------------------------------------------------
c LSD potential from [c](A1)
c ECRS = dEc/drs [c](A2)
c ECZET=dEc/dzeta [c](A3)
c FZ = dF/dzeta [c](A4)
      ECRS = EURS*(1.D0-F*Z4)+EPRS*F*Z4-ALFRSM*F*(1.D0-Z4)/FZZ
      FZ = THRD4*((1.D0+ZET)**THRD-(1.D0-ZET)**THRD)/GAM
      ECZET = 4.D0*(ZET**3)*F*(EP-EU+ALFM/FZZ)+FZ*(Z4*EP-Z4*EU
     1        -(1.D0-Z4)*ALFM/FZZ)
      COMM = EC -RS*ECRS/3.D0-ZET*ECZET
      VCUP = COMM + ECZET
      VCDN = COMM - ECZET
      if(lgga.eq.0)return
c----------------------------------------------------------------------
c----------------------------------------------------------------------
c PBE correlation energy
c G=phi(zeta), given after [a](3)
c DELT=bet/gamma
c B=A of [a](8)
      G=((1.d0+ZET)**thrd2+(1.d0-ZET)**thrd2)/2.d0
      G3 = G**3
      PON=-EC/(G3*gamma)
      B = DELT/(DEXP(PON)-1.D0)
      B2 = B*B
      T2 = T*T
      T4 = T2*T2
      RS2 = RS*RS
      RS3 = RS2*RS
      Q4 = 1.D0+B*T2
      Q5 = 1.D0+B*T2+B2*T4
      H = G3*(BET/DELT)*DLOG(1.D0+DELT*Q4*T2/Q5)
      if(lpot.eq.0)return
c----------------------------------------------------------------------
c----------------------------------------------------------------------
C ENERGY DONE. NOW THE POTENTIAL, using appendix E of [b].
      G4 = G3*G
      T6 = T4*T2
      RSTHRD = RS/3.D0
      GZ=(((1.d0+zet)**2+eta)**sixthm-
     1((1.d0-zet)**2+eta)**sixthm)/3.d0
      FAC = DELT/B+1.D0
      BG = -3.D0*B2*EC*FAC/(BET*G4)
      BEC = B2*FAC/(BET*G3)
      Q8 = Q5*Q5+DELT*Q4*Q5*T2
      Q9 = 1.D0+2.D0*B*T2
      hB = -BET*G3*B*T6*(2.D0+B*T2)/Q8
      hRS = -RSTHRD*hB*BEC*ECRS
      FACT0 = 2.D0*DELT-6.D0*B
      FACT1 = Q5*Q9+Q4*Q9*Q9
      hBT = 2.D0*BET*G3*T4*((Q4*Q5*FACT0-DELT*FACT1)/Q8)/Q8
      hRST = RSTHRD*T2*hBT*BEC*ECRS
      hZ = 3.D0*GZ*h/G + hB*(BG*GZ+BEC*ECZET)
      hT = 2.d0*BET*G3*Q9/Q8
      hZT = 3.D0*GZ*hT/G+hBT*(BG*GZ+BEC*ECZET)
      FACT2 = Q4*Q5+B*T2*(Q4*Q9+Q5)
      FACT3 = 2.D0*B*Q5*Q9+DELT*FACT2
      hTT = 4.D0*BET*G3*T*(2.D0*B/Q8-(Q9*FACT3/Q8)/Q8)
      COMM = H+HRS+HRST+T2*HT/6.D0+7.D0*T2*T*HTT/6.D0
      PREF = HZ-GZ*T2*HT/G
      FACT5 = GZ*(2.D0*HT+T*HTT)/G
      COMM = COMM-PREF*ZET-UU*HTT-VV*HT-WW*(HZT-FACT5)
      DVCUP = COMM + PREF
      DVCDN = COMM - PREF
      RETURN
      END
c----------------------------------------------------------------------
c######################################################################
c----------------------------------------------------------------------
      SUBROUTINE GCOR2(A,A1,B1,B2,B3,B4,rtrs,GG,GGRS)
c slimmed down version of GCOR used in PW91 routines, to interpolate
c LSD correlation energy, as given by (10) of
c J. P. Perdew and Y. Wang, Phys. Rev. B {\bf 45}, 13244 (1992).
c K. Burke, May 11, 1996.
      IMPLICIT REAL*8 (A-H,O-Z)
      Q0 = -2.D0*A*(1.D0+A1*rtrs*rtrs)
      Q1 = 2.D0*A*rtrs*(B1+rtrs*(B2+rtrs*(B3+B4*rtrs)))
      Q2 = DLOG(1.D0+1.D0/Q1)
      GG = Q0*Q2
      Q3 = A*(B1/rtrs+2.D0*B2+rtrs*(3.D0*B3+4.D0*B4*rtrs))
      GGRS = -2.D0*A*A1*Q2-Q0*Q3/(Q1*(1.d0+Q1))
      RETURN
      END
c----------------------------------------------------------------------
c######################################################################
c----------------------------------------------------------------------
      SUBROUTINE EXCHPW91(D,S,U,V,EX,VX)
C  GGA91 EXCHANGE FOR A SPIN-UNPOLARIZED ELECTRONIC SYSTEM
C  INPUT D : DENSITY
C  INPUT S:  ABS(GRAD D)/(2*KF*D)
C  INPUT U:  (GRAD D)*GRAD(ABS(GRAD D))/(D**2 * (2*KF)**3)
C  INPUT V: (LAPLACIAN D)/(D*(2*KF)**2)
C  OUTPUT:  EXCHANGE ENERGY PER ELECTRON (EX) AND POTENTIAL (VX)
      IMPLICIT REAL*8 (A-H,O-Z)
      parameter(a1=0.19645D0,a2=0.27430D0,a3=0.15084D0,a4=100.d0)
      parameter(ax=-0.7385588D0,a=7.7956D0,b1=0.004d0)
      parameter(thrd=0.333333333333D0,thrd4=1.33333333333D0)
c for Becke exchange, set a3=b1=0
      FAC = AX*D**THRD
      S2 = S*S
      S3 = S2*S
      S4 = S2*S2
      P0 = 1.D0/DSQRT(1.D0+A*A*S2)
      P1 = DLOG(A*S+1.D0/P0)
      P2 = DEXP(-A4*S2)
      P3 = 1.D0/(1.D0+A1*S*P1+B1*S4)
      P4 = 1.D0+A1*S*P1+(A2-A3*P2)*S2
      F = P3*P4
      EX = FAC*F
C  LOCAL EXCHANGE OPTION
C     EX = FAC
C  ENERGY DONE. NOW THE POTENTIAL:
      P5 = B1*S2-(A2-A3*P2)
      P6 = A1*S*(P1+A*S*P0)
      P7 = 2.D0*(A2-A3*P2)+2.D0*A3*A4*S2*P2-4.D0*B1*S2*F
      FS = P3*(P3*P5*P6+P7)
      P8 = 2.D0*S*(B1-A3*A4*P2)
      P9 = A1*P1+A*A1*S*P0*(3.D0-A*A*S2*P0*P0)
      P10 = 4.D0*A3*A4*S*P2*(2.D0-A4*S2)-8.D0*B1*S*F-4.D0*B1*S3*FS
      P11 = -P3*P3*(A1*P1+A*A1*S*P0+4.D0*B1*S3)
      FSS = P3*P3*(P5*P9+P6*P8)+2.D0*P3*P5*P6*P11+P3*P10+P7*P11
      VX = FAC*(THRD4*F-(U-THRD4*S3)*FSS-V*FS)
C  LOCAL EXCHANGE OPTION:
C     VX = FAC*THRD4
      RETURN
      END
c----------------------------------------------------------------------
c######################################################################
c----------------------------------------------------------------------
      SUBROUTINE CORLSD(RS,ZET,EC,VCUP,VCDN,ECRS,ECZET,ALFC)
C  UNIFORM-GAS CORRELATION OF PERDEW AND WANG 1991
C  INPUT: SEITZ RADIUS (RS), RELATIVE SPIN POLARIZATION (ZET)
C  OUTPUT: CORRELATION ENERGY PER ELECTRON (EC), UP- AND DOWN-SPIN
C     POTENTIALS (VCUP,VCDN), DERIVATIVES OF EC WRT RS (ECRS) & ZET (ECZET)
C  OUTPUT: CORRELATION CONTRIBUTION (ALFC) TO THE SPIN STIFFNESS
      IMPLICIT REAL*8 (A-H,O-Z)
      parameter(gam=0.5198421D0,fzz=1.709921D0)
      parameter(thrd=0.333333333333D0,thrd4=1.333333333333D0)
      F = ((1.D0+ZET)**THRD4+(1.D0-ZET)**THRD4-2.D0)/GAM
      CALL GCOR(0.0310907D0,0.21370D0,7.5957D0,3.5876D0,1.6382D0,
     1    0.49294D0,1.00D0,RS,EU,EURS)
      CALL GCOR(0.01554535D0,0.20548D0,14.1189D0,6.1977D0,3.3662D0,
     1    0.62517D0,1.00D0,RS,EP,EPRS)
      CALL GCOR(0.0168869D0,0.11125D0,10.357D0,3.6231D0,0.88026D0,
     1    0.49671D0,1.00D0,RS,ALFM,ALFRSM)
C  ALFM IS MINUS THE SPIN STIFFNESS ALFC
      ALFC = -ALFM
      Z4 = ZET**4
      EC = EU*(1.D0-F*Z4)+EP*F*Z4-ALFM*F*(1.D0-Z4)/FZZ
C  ENERGY DONE. NOW THE POTENTIAL:
      ECRS = EURS*(1.D0-F*Z4)+EPRS*F*Z4-ALFRSM*F*(1.D0-Z4)/FZZ
      FZ = THRD4*((1.D0+ZET)**THRD-(1.D0-ZET)**THRD)/GAM
      ECZET = 4.D0*(ZET**3)*F*(EP-EU+ALFM/FZZ)+FZ*(Z4*EP-Z4*EU
     1        -(1.D0-Z4)*ALFM/FZZ)
      COMM = EC -RS*ECRS/3.D0-ZET*ECZET
      VCUP = COMM + ECZET
      VCDN = COMM - ECZET
      RETURN
      END
c----------------------------------------------------------------------
c######################################################################
c----------------------------------------------------------------------
      SUBROUTINE GCOR(A,A1,B1,B2,B3,B4,P,RS,GG,GGRS)
C  CALLED BY SUBROUTINE CORLSD
      IMPLICIT REAL*8 (A-H,O-Z)
      P1 = P + 1.D0
      Q0 = -2.D0*A*(1.D0+A1*RS)
      RS12 = DSQRT(RS)
      RS32 = RS12**3
      RSP = RS**P
      Q1 = 2.D0*A*(B1*RS12+B2*RS+B3*RS32+B4*RS*RSP)
      Q2 = DLOG(1.D0+1.D0/Q1)
      GG = Q0*Q2
      Q3 = A*(B1/RS12+2.D0*B2+3.D0*B3*RS12+2.D0*B4*P1*RSP)
      GGRS = -2.D0*A*A1*Q2-Q0*Q3/(Q1**2+Q1)
      RETURN
      END
c----------------------------------------------------------------------
c######################################################################
c----------------------------------------------------------------------
      SUBROUTINE CORpw91(RS,ZET,G,EC,ECRS,ECZET,T,UU,VV,WW,H,
     1                   DVCUP,DVCDN)
C  pw91 CORRELATION, modified by K. Burke to put all arguments
c  as variables in calling statement, rather than in common block
c  May, 1996.
C  INPUT RS: SEITZ RADIUS
C  INPUT ZET: RELATIVE SPIN POLARIZATION
C  INPUT T: ABS(GRAD D)/(D*2.*KS*G)
C  INPUT UU: (GRAD D)*GRAD(ABS(GRAD D))/(D**2 * (2*KS*G)**3)
C  INPUT VV: (LAPLACIAN D)/(D * (2*KS*G)**2)
C  INPUT WW:  (GRAD D)*(GRAD ZET)/(D * (2*KS*G)**2
C  OUTPUT H: NONLOCAL PART OF CORRELATION ENERGY PER ELECTRON
C  OUTPUT DVCUP,DVCDN:  NONLOCAL PARTS OF CORRELATION POTENTIALS
      IMPLICIT REAL*8 (A-H,O-Z)
      parameter(xnu=15.75592D0,cc0=0.004235D0,cx=-0.001667212D0)
      parameter(alf=0.09D0)
      parameter(c1=0.002568D0,c2=0.023266D0,c3=7.389D-6,c4=8.723D0)
      parameter(c5=0.472D0,c6=7.389D-2,a4=100.D0)
      parameter(thrdm=-0.333333333333D0,thrd2=0.666666666667D0)
      BET = XNU*CC0
      DELT = 2.D0*ALF/BET
      G3 = G**3
      G4 = G3*G
      PON = -DELT*EC/(G3*BET)
      B = DELT/(DEXP(PON)-1.D0)
      B2 = B*B
      T2 = T*T
      T4 = T2*T2
      T6 = T4*T2
      RS2 = RS*RS
      RS3 = RS2*RS
      Q4 = 1.D0+B*T2
      Q5 = 1.D0+B*T2+B2*T4
      Q6 = C1+C2*RS+C3*RS2
      Q7 = 1.D0+C4*RS+C5*RS2+C6*RS3
      CC = -CX + Q6/Q7
      R0 = 0.663436444d0*rs
      R1 = A4*R0*G4
      COEFF = CC-CC0-3.D0*CX/7.D0
      R2 = XNU*COEFF*G3
      R3 = DEXP(-R1*T2)
      H0 = G3*(BET/DELT)*DLOG(1.D0+DELT*Q4*T2/Q5)
      H1 = R3*R2*T2
      H = H0+H1
C  LOCAL CORRELATION OPTION:
C     H = 0.0D0
C  ENERGY DONE. NOW THE POTENTIAL:
      CCRS = (C2+2.*C3*RS)/Q7 - Q6*(C4+2.*C5*RS+3.*C6*RS2)/Q7**2
      RSTHRD = RS/3.D0
      R4 = RSTHRD*CCRS/COEFF
      GZ = ((1.D0+ZET)**THRDM - (1.D0-ZET)**THRDM)/3.D0
      FAC = DELT/B+1.D0
      BG = -3.D0*B2*EC*FAC/(BET*G4)
      BEC = B2*FAC/(BET*G3)
      Q8 = Q5*Q5+DELT*Q4*Q5*T2
      Q9 = 1.D0+2.D0*B*T2
      H0B = -BET*G3*B*T6*(2.D0+B*T2)/Q8
      H0RS = -RSTHRD*H0B*BEC*ECRS
      FACT0 = 2.D0*DELT-6.D0*B
      FACT1 = Q5*Q9+Q4*Q9*Q9
      H0BT = 2.D0*BET*G3*T4*((Q4*Q5*FACT0-DELT*FACT1)/Q8)/Q8
      H0RST = RSTHRD*T2*H0BT*BEC*ECRS
      H0Z = 3.D0*GZ*H0/G + H0B*(BG*GZ+BEC*ECZET)
      H0T = 2.*BET*G3*Q9/Q8
      H0ZT = 3.D0*GZ*H0T/G+H0BT*(BG*GZ+BEC*ECZET)
      FACT2 = Q4*Q5+B*T2*(Q4*Q9+Q5)
      FACT3 = 2.D0*B*Q5*Q9+DELT*FACT2
      H0TT = 4.D0*BET*G3*T*(2.D0*B/Q8-(Q9*FACT3/Q8)/Q8)
      H1RS = R3*R2*T2*(-R4+R1*T2/3.D0)
      FACT4 = 2.D0-R1*T2
      H1RST = R3*R2*T2*(2.D0*R4*(1.D0-R1*T2)-THRD2*R1*T2*FACT4)
      H1Z = GZ*R3*R2*T2*(3.D0-4.D0*R1*T2)/G
      H1T = 2.D0*R3*R2*(1.D0-R1*T2)
      H1ZT = 2.D0*GZ*R3*R2*(3.D0-11.D0*R1*T2+4.D0*R1*R1*T4)/G
      H1TT = 4.D0*R3*R2*R1*T*(-2.D0+R1*T2)
      HRS = H0RS+H1RS
      HRST = H0RST+H1RST
      HT = H0T+H1T
      HTT = H0TT+H1TT
      HZ = H0Z+H1Z
      HZT = H0ZT+H1ZT
      COMM = H+HRS+HRST+T2*HT/6.D0+7.D0*T2*T*HTT/6.D0
      PREF = HZ-GZ*T2*HT/G
      FACT5 = GZ*(2.D0*HT+T*HTT)/G
      COMM = COMM-PREF*ZET-UU*HTT-VV*HT-WW*(HZT-FACT5)
      DVCUP = COMM + PREF
      DVCDN = COMM - PREF
C  LOCAL CORRELATION OPTION:
C     DVCUP = 0.0D0
C     DVCDN = 0.0D0
      RETURN
      END

max_stars_repo_name: qsnake/petot

max_issues_repo_name: qsnake/petot

max_forks_repo_name: qsnake/petot


text:
c
c     Matrix operation library
c
      subroutine v_clear(a,n)
      implicit real*8(a-h,o-z)
      real*8 a(n)
      a=0.d0
c      do 10 i=1,n
c 10      a(i)=0.d0
      return
      end
c ------------------------------------------------------------------
      subroutine m_clear(a,n)
      implicit real*8(a-h,o-z)
      real*8 a(*)
      a(1:n*n)=0.d0
c      do 10 i=1,n*n
c 10      a(i)=0.d0
      return
      end
c ------------------------------------------------------------------
      subroutine m_unit(a,n)
      implicit real*8(a-h,o-z)
      real*8 a(n,n)
      a=0.d0
c      do 20 i=1,n
c         do 20 j=1,n
c            a(j,i)=0.
c 20   continue
      do i=1,n
        a(i,i)=1.d0
      enddo
      return
      end
c ------------------------------------------------------------------
      subroutine v_copy(a,b,n)
      implicit real*8(a-h,o-z)
      real*8 a(n),b(n)
      b=a
c      do 10 i=1,n
c 10      b(i)=a(i)
      return
      end
c ------------------------------------------------------------------
      subroutine m_copy_6d(a,b)
      real*8 a(36),b(36)
      b=a
c      do 10 i=1,36
c 10      b(i)=a(i)
      return
      end
c ------------------------------------------------------------------
      subroutine v_chg_sign(a,n)
      real*8 a(n)
      a=-a
c      do 10 i=1,n
c 10      a(i)=-a(i)
      return
      end
c ------------------------------------------------------------------
      subroutine m_print(a,n)
      integer*4 ,intent(in):: n
      real*8 ,intent(in):: a(n,n)
      integer*4 i
c      write(*,*) 'Matrix print out : Ndim=',n
      do i=1,n
        write(*,'(1P,(6D12.5))') (a(i,j),j=1,n)
      enddo
      return
      end
c ------------------------------------------------------------------
      subroutine m_add_6d(a,b,c)
      implicit real*8(a-h,o-z)
      real*8 a(36),b(36),c(36)
      c=a+b
c      do 10 i=1,36
c      c(i)=a(i)+b(i)
c 10   continue
      return
      end
c
c ------------------------------------------------------------------
      subroutine m_sub_6d(a,b,c)
      implicit real*8(a-h,o-z)
      real*8 a(36),b(36),c(36)
      c=a-b
c      do 10 i=1,36
c      c(i)=a(i)-b(i)
c 10   continue
      return
      end
c
c ------------------------------------------------------------------
      subroutine m_mul_6d(a,b,c)
      implicit real*8(a-h,o-z)
      real*8 a(6,6),b(6,6),c(6,6)
c      real*8 t(6,6)
c
      c=matmul(a,b)
c      call m_clear(t,6)
c      do 10 i=1,6
c      do 20 j=1,6
c      t(i,j)=0.d0
c      do 30 k=1,6
c         t(i,j)=t(i,j)+a(i,k)*b(k,j)
c 30   continue
c 20   continue
c 10   continue
c
c     do 40 i=1,6
c      do 40 j=1,6
c      c(i,j)=t(i,j)
c 40   continue
      return
      end
c
c ------------------------------------------------------------------
      subroutine m_mul3_6d(a,b,c,d)
      implicit real*8(a-h,o-z)
      real*8 a(*),b(*),c(*),d(*)
      real*8 t(6,6)
c
      call m_mul_6d(a,b,t)
      call m_mul_6d(t,c,d)
      return
      end
c
c ------------------------------------------------------------------
      subroutine m_transpose_6d(a,b)
      implicit real*8(a-h,o-z)
      real*8 a(6,6),b(6,6)
c      real*8 t(6,6)
      b=transpose(a)
c
c      do 10 i=1,6
c      do 10 j=1,6
c 10      t(i,j)=a(j,i)
c      do 20 i=1,6
c      do 20 j=1,6
c 20      b(j,i)=t(j,i)
c      return
      end
c ------------------------------------------------------------------
      subroutine m_sim_tr_6d(a,b,c)
      implicit real*8(a-h,o-z)
      real*8 a(*),b(*),c(*)
      real*8 t(6,6)
c
      call m_transpose_6d(a,t)
      call m_mul3_6d(a,b,t,c)
      return
      end
c ------------------------------------------------------------------
      subroutine m_mul_5d(a,b,c)
      implicit real*8(a-h,o-z)
      real*8 a(5,5),b(5,5),c(5,5)
c      real*8 t(5,5)
c
      c=matmul(a,b)
c$$$      call m_clear(t,5)
c$$$      do 10 i=1,5
c$$$      do 20 j=1,5
c$$$      t(i,j)=0.d0
c$$$      do 30 k=1,5
c$$$         t(i,j)=t(i,j)+a(i,k)*b(k,j)
c$$$ 30   continue
c$$$ 20   continue
c$$$ 10   continue
c$$$c
c$$$      do 40 i=1,5
c$$$      do 40 j=1,5
c$$$      c(i,j)=t(i,j)
c$$$ 40   continue
      return
      end
c
c ------------------------------------------------------------------
      subroutine m_mul3_5d(a,b,c,d)
      implicit real*8(a-h,o-z)
      real*8 a(*),b(*),c(*),d(*)
      real*8 t(5,5)
c
      call m_mul_5d(a,b,t)
      call m_mul_5d(t,c,d)
      return
      end
c
c ------------------------------------------------------------------
      subroutine m_transpose_5d(a,b)
      implicit real*8(a-h,o-z)
      real*8 a(5,5),b(5,5)
c      real*8 t(5,5)
c
      b=transpose(a)
c$$$      do 10 i=1,5
c$$$      do 10 j=1,5
c$$$ 10      t(i,j)=a(j,i)
c$$$      do 20 i=1,5
c$$$      do 20 j=1,5
c$$$ 20      b(j,i)=t(j,i)
      return
      end
c ------------------------------------------------------------------
      subroutine m_sim_tr_5d(a,b,c)
      implicit real*8(a-h,o-z)
      real*8 a(*),b(*),c(*)
      real*8 t(6,6)
c
      call m_transpose_5d(a,t)
      call m_mul3_5d(a,b,t,c)
      return
      end
c ------------------------------------------------------------------
      subroutine m_symp_set_6d(S)
      real*8 S(6,6)
c
c      call m_clear(S,6)
      s=0.d0
      S(1,2)=1.d0
      S(2,1)=-1.d0
      S(3,4)=1.d0
      S(4,3)=-1.d0
      S(5,6)=1.d0
      S(6,5)=-1.d0
      return
      end
c-----------------------------------------------------------------
      subroutine m_symmet_set_6d(p_in,blist)
      real*8 p_in(*),blist(*)
      k=1
      do i=1,6
       do j=1,6
         if(i.le.j) then
            blist(i+(j-1)*6)=p_in(k)
            k=k+1
         else
            blist(i+(j-1)*6)=blist(j+(i-1)*6)
         endif
       enddo
      enddo

      return
      end
c ------------------------------------------------------------------
      subroutine is_symplectic_6d(a)
      real*8 a(6,6),t(6,6),S(6,6)
      call m_symp_set_6d(S)
      call m_transpose_6d(a,t)
      call m_mul3_6d(t,S,a,t)
      write(*,*) ' Symplectic check'
      call m_print(t,6)
      return
      end
c ------------------------------------------------------------------
      subroutine m_B_set_6d(ax,bx,ay,by,az,bz,m_B)
      implicit real*8(a-h,o-z)
      real*8 ax,bx,ay,by,az,bz,m_B(6,6)
c
      call m_clear(m_B,6)
      sqrb=sqrt(bx)
      m_B(1,1)=sqrb
      m_B(2,1)=-ax/sqrb
      m_B(2,2)=1.d0/sqrb
      sqrb=sqrt(by)
      m_B(3,3)=sqrb
      m_B(4,3)=-ay/sqrb
      m_B(4,4)=1.d0/sqrb
      sqrb=sqrt(bz)
      m_B(5,5)=sqrb
      m_B(6,5)=-az/sqrb
      m_B(6,6)=1.d0/sqrb
      return
      end
c ------------------------------------------------------------------
      subroutine m_R_set_6d(r1,r2,r3,r4,m_R)
      implicit real*8(a-h,o-z)
      real*8 r1,r2,r3,r4,m_R(6,6)
      real*8 rmu
      rmu=sqrt(1.d0-r1*r4+r2*r3)
c
      call m_clear(m_R,6)
      m_R(1,1)=rmu
      m_R(2,2)=rmu
      m_R(3,3)=rmu
      m_R(4,4)=rmu
      m_R(5,5)=1.d0
      m_R(6,6)=1.d0
c
      m_R(3,1)=-r1
      m_R(3,2)=-r2
      m_R(4,1)=-r3
      m_R(4,2)=-r4
      m_R(1,3)=r4
      m_R(1,4)=-r2
      m_R(2,3)=-r3
      m_R(2,4)=r1
c
      return
      end
c ------------------------------------------------------------------
      subroutine m_H_set_6d(ex,epx,ey,epy,zx,zpx,zy,zpy,m_H)
      implicit real*8(a-h,o-z)
      real*8 ex,epx,ey,epy,zx,zpx,zy,zpy,m_H(6,6)
c
      call m_clear(m_H,6)
      detHx=zx*epx-zpx*ex
      detHy=zy*epy-zpy*ey
      a2=1.d0-detHx-detHy
      if(a2.gt.0.) then
        a=sqrt(a2)
      else
        write(*,*) ' det Hx +det Hy >1 ',detHx,detHy
        stop
      endif
      ai=1.d0/(1.d0+a)
c
      m_H(1,1)=1.d0-detHx*ai
      m_H(2,2)=m_H(1,1)
      m_H(3,3)=1.d0-detHy*ai
      m_H(4,4)=m_H(3,3)
      m_H(5,5)=a
      m_H(6,6)=m_H(5,5)
c
      m_H(1,5)=zx
      m_H(1,6)=ex
      m_H(2,5)=zpx
      m_H(2,6)=epx
      m_H(3,5)=zy
      m_H(3,6)=ey
      m_H(4,5)=zpy
      m_H(4,6)=epy
c
      m_H(5,1)=-epx
      m_H(5,2)=ex
      m_H(6,1)=zpx
      m_H(6,2)=-zx
      m_H(5,3)=-epy
      m_H(5,4)=ey
      m_H(6,3)=zpy
      m_H(6,4)=-zy
c
      m_H(1,3)=(ex*zpy-epy*zx)*ai
      m_H(1,4)=(ey*zx-ex*zy)*ai
      m_H(2,3)=(epx*zpy-epy*zpx)*ai
      m_H(2,4)=(ey*zpx-epx*zy)*ai
      m_H(3,1)=(ey*zpx-epx*zy)*ai
      m_H(3,2)=(-ey*zx+ex*zy)*ai
      m_H(4,1)=(epy*zpx-epx*zpy)*ai
      m_H(4,2)=(ex*zpy-epy*zx)*ai
      return
      end
c ------------------------------------------------------------------
      subroutine m_NtoP_set_6d(a,m_NtoP)
      implicit real*8(a-h,o-z)
      real*8 a(*),m_NtoP(6,6)
      real*8 t(6,6)
c
c      call m_clear(m_NtoP,6)
       if(a(2).le.0.) then
          write(*,*) ' Finite BETAX is required.'
          write(*,*) ' BETAX is set to be 1..'
          a(2)=1.d0
       endif
       if(a(5).le.0.) then
          write(*,*) ' Finite BETAY is required.'
          write(*,*) ' BETAY is set to be 1..'
          a(5)=1.d0
       endif
       if(a(27).le.0.) then
          write(*,*) ' Finite BETAZ is required.'
          write(*,*) ' BETAZ is set to be 1..'
          a(27)=1.d0
       endif
c
       call m_B_set_6d(a(1),a(2),a(4),a(5),a(22),a(27),m_NtoP)
       call m_R_set_6d(a(7),a(8),a(9),a(10),t)
       call m_mul_6d(t,m_NtoP,m_NtoP)
       call m_H_set_6d(a(11),a(12),a(13),a(14),
     +      a(15),a(16),a(17),a(18),t)
       call m_mul_6d(t,m_NtoP,m_NtoP)
c       write(*,*) ' m_NtoP_set_6d'
c       call m_print(m_NtoP,6)
      return
      end
c ------------------------------------------------------------------
      subroutine m_PtoN_set_6d(a,m_PtoN)
      implicit real*8(a-h,o-z)
      real*8 a(*),m_PtoN(6,6)
c
      call m_NtoP_set_6d(a,m_PtoN)
      call m_sym_inv_6d(m_PtoN,m_PtoN)
      return
      end
c ------------------------------------------------------------------
      subroutine m_sym_inv_6d(a,b)
      implicit real*8(a-h,o-z)
      real*8 a(6,6),b(6,6)
      real*8 S6(6,6)
c
      call m_symp_set_6d(S6)
      call m_sim_tr_6d(S6,a,b)
      return
      end
c ------------------------------------------------------------------
      subroutine m_prot_set_6d(p,a)
      use macmath
      real*8 px,py,pz,p(*),a(6,6)
      call m_clear(a,6)
      px=p(3)*2.d0*pi
      py=p(6)*2.d0*pi
      pz=p(25)*2.d0*pi
      a(1,1)=cos(px)
      a(1,2)=sin(px)
      a(2,1)=-sin(px)
      a(2,2)=cos(px)
      a(3,3)=cos(py)
      a(3,4)=sin(py)
      a(4,3)=-sin(py)
      a(4,4)=cos(py)
      a(5,5)=cos(pz)
      a(5,6)=-sin(pz)
      a(6,5)=sin(pz)
      a(6,6)=cos(pz)
      return
      end
c ------------------------------------------------------------------
      subroutine m_Emit_set_6d(emitx,emity,emitz,m_Emit)
      implicit real*8(a-h,o-z)
      real*8 emitx,emity,emitz,m_Emit(6,6)
c
      call m_clear(m_Emit,6)
      m_Emit(1,1)=emitx
      m_Emit(2,2)=emitx
      m_Emit(3,3)=emity
      m_Emit(4,4)=emity
      m_Emit(5,5)=emitz
      m_Emit(6,6)=emitz
      return
      end
c ------------------------------------------------------------------
      subroutine m_Crs_set_6d(x_angle,m_Crs)
      implicit real*8(a-h,o-z)
      real*8 x_angle,m_Crs(6,6)
c
      call m_clear(m_Crs,6)
      cosxi=1.d0/cos(x_angle)
      tanx=tan(x_angle)
      m_Crs(1,1)=1.d0
      m_Crs(2,2)=cosxi
      m_Crs(3,3)=1.d0
      m_Crs(4,4)=cosxi
      m_Crs(5,5)=cosxi
      m_Crs(6,6)=1.d0
c
      m_Crs(1,5)=tanx
      m_Crs(6,2)=-tanx
      return
      end
c ------------------------------------------------------------------
      subroutine m_inverse_6d(a,b)
      real*8 a(6,6),b(6,6),t(6,6)
      integer indx(6)
c
      t=a
c      do 10 i=1,6
c      do 10 j=1,6
c 10     t(j,i)=a(j,i)
      call m_inverse(t,6,6,indx,b)
      return
      end
c ------------------------------------------------------------------
      subroutine m_zslice_5d(x,a,b,c)
      real*8 x(6,6),a(5,5),b(5),c
c
      a(1:4,1:4)=x(1:4,1:4)
c      do 10 j=1,4
c      do 10 i=1,4
c 10      a(i,j)=x(i,j)
      do i=1,4
        b(i)=x(5,i)
        a(5,i)=x(6,i)
        a(i,5)=x(i,6)
      enddo
      c=x(5,5)
      b(5)=x(5,6)
      a(5,5)=x(6,6)
      return
      end
c ------------------------------------------------------------------
      subroutine m_normalize_5d(a,b)
      real*8 a(5,5),b(10)
      real*8 ta(5,5),tb(10),rnorm,remax,absa
      integer naxis(5),ikeep(5)
      ikeep=0
c      do 5 i=1,5
c 5       ikeep(i)=0
c
      do i=1,5
        rnorm=0.d0
        remax=0.
        naxis(i)=0
        do j=1,5
          rnorm=rnorm+a(j,i)**2
        enddo
        rnorm=1.d0/sqrt(rnorm)
        do j=1,5
          ta(j,i)=a(j,i)*rnorm
          absa=abs(ta(j,i))
          if(absa.gt.remax.and.ikeep(j).eq.0) then
            remax=absa
            naxis(i)=j
          endif
        enddo
        ikeep(naxis(i))=1
      enddo
      do i=1,5
        if(ikeep(i).eq.0) then
          write(*,*) 'matrix.f : Normal axis is determined'
          write(*,*) naxis
          stop
        endif
      enddo
      do i=1,5
        k=naxis(i)
        do j=1,5
          a(j,k)=ta(j,i)
        enddo
        tb(2*k-1)=b(2*i-1)
        tb(2*k)=b(2*i)
      enddo
      b=tb
c      do 60 i=1,10
c 60     b(i)=tb(i)
      return
      end
c ------------------------------------------------------------------
      subroutine m_normalize_6d(a,b)
      real*8 a(6,6),b(12)
      real*8 ta(6,6),tb(12),rnorm,remax,absa
      integer naxis(6),ikeep(6)
      ikeep=0
c      do 5 i=1,6
c 5       ikeep(i)=0
c
      do i=1,6
        rnorm=0.d0
        remax=0.
        naxis(i)=0
        do j=1,6
          rnorm=rnorm+a(j,i)**2
        enddo
        rnorm=1.d0/sqrt(rnorm)
        do j=1,6
          ta(j,i)=a(j,i)*rnorm
          absa=abs(ta(j,i))
          if(absa.gt.remax.and.ikeep(j).eq.0) then
            remax=absa
            naxis(i)=j
          endif
        enddo
        ikeep(naxis(i))=1
      enddo
      do i=1,6
        if(ikeep(i).eq.0) then
          write(*,*) 'matrix.f : Normal axis is determined'
          write(*,*) naxis
          stop
        endif
      enddo
      do i=1,6
        k=naxis(i)
        do j=1,6
          a(j,k)=ta(j,i)
        enddo
        tb(2*k-1)=b(2*i-1)
        tb(2*k)=b(2*i)
      enddo
      b=tb
c      do 60 i=1,12
c 60     b(i)=tb(i)
        return
        end
c ------------------------------------------------------------------
      subroutine m_diag_5d(a1,a2,a3,a4,a5,b)
      real*8 a1,a2,a3,a4,a5,b(5,5)
      call m_clear(b,5)
      b(1,1)=a1
      b(2,2)=a2
      b(3,3)=a3
      b(4,4)=a4
      b(5,5)=a5
      return
      end
c ------------------------------------------------------------------
      subroutine mv_mul_5d(a,b,c)
      real*8 a(5,5),b(5),c(5)
      c=matmul(a,b)
c      do 10 i=1,5
c       t(i)=0.d0
c      do 10 j=1,5
c 10    t(i)=t(i)+a(i,j)*b(j)
c      do 20 i=1,5
c 20      c(i)=t(i)
      return
      end
c ------------------------------------------------------------------
      subroutine mv_mul_6d(a,b,c)
      real*8 a(6,6),b(6),c(6)
      c=matmul(a,b)
c$$$      do 10 i=1,6
c$$$       t(i)=0.d0
c$$$      do 10 j=1,6
c$$$ 10    t(i)=t(i)+a(i,j)*b(j)
c$$$      do 20 i=1,6
c$$$ 20      c(i)=t(i)
      return
      end
c ------------------------------------------------------------------
      subroutine benv_norm_axis(env,vec,t_angle)
      implicit real*8(a-h,o-z)
      real*8 env(5,5),vec(5),t_angle
      real*8 u(5,5),tmp(5,5),xx(5)
      t_angle=atan2(2.d0*env(1,3),env(1,1)-env(3,3))*0.5d0
      cost=cos(t_angle)
      sint=sin(t_angle)
      call m_clear(u,5)
      u(1,1)=cost
      u(2,2)=cost
      u(3,3)=cost
      u(4,4)=cost
      u(5,5)=1.d0
      u(1,3)=sint
      u(2,4)=sint
      u(3,1)=-sint
      u(4,2)=-sint
c
      call m_clear(tmp,5)
      do i=1,5
        do j=1,5
          do k=1,5
            tmp(i,j)=tmp(i,j)+env(i,k)*u(j,k)
          enddo
        enddo
      enddo
      call m_clear(env,5)
      do i=1,5
        do j=1,5
          do k=1,5
            env(i,j)=env(i,j)+u(i,k)*tmp(k,j)
          enddo
        enddo
      enddo
      do i=1,5
        xx(i)=0.
        do j=1,5
          xx(i)=xx(i)+u(i,j)*vec(j)
        enddo
      enddo
      vec=xx
c      do 40 i=1,5
c 40      vec(i)=xx(i)
      return
      end
c ------------------------------------------------------------------
      subroutine m_inverse(a,n,np,indx,y)
      integer n,np,indx(np)
      real*8 a(np,np),y(np,np),d
      call m_unit(y,np)
      call ludcmp(a,n,np,indx,d)
      do 10 j=1,n
      call lubksb(a,n,np,indx,y(1,j))
  10  continue
      return
      end
c ------------------------------------------------------------------
      subroutine ludcmp(a,n,np,indx,d)
      implicit real*8(a-h,o-z)
      parameter (nmax=100,tiny=1.e-20)
      integer n,np,indx(np)
      real*8 a(np,np),d,vv(nmax)
c     begin initialize for preventing compiler warning
      imax=0
c     end   initialize for preventing compiler warning
      d=1.d0
      do 10 i=1,n
        aamax=0.d0
        do 11 j=1,n
          if(abs(a(i,j)).gt.aamax) aamax=abs(a(i,j))
  11    continue
        if(aamax.eq.0.d0) WRITE(*,*) ' Singular matrix '
        vv(i)=1./aamax
  10  continue
      do 19 j=1,n
        if(j.gt.1) then
          do 14 i=1,j-1
            sum=a(i,j)
            if(i.gt.1) then
              do 13 k=1,i-1
                sum=sum-a(i,k)*a(k,j)
  13          continue
              a(i,j)=sum
            endif
  14      continue
        endif
        aamax=0.d0
        do 16 i=j,n
          sum=a(i,j)
          if(j.gt.1) then
            do 15 k=1,j-1
              sum=sum-a(i,k)*a(k,j)
  15        continue
            a(i,j)=sum
          endif
          dum=vv(i)*abs(sum)
          if(dum.gt.aamax) then
            imax=i
            aamax=dum
          endif
  16    continue
        if(j.ne.imax) then
          do 17 k=1,n
            dum=a(imax,k)
            a(imax,k)=a(j,k)
            a(j,k)=dum
  17      continue
          d=-d
          vv(imax)=vv(j)
        endif
        indx(j)=imax
        if(j.ne.n) then
          if(a(j,j).eq.0.d0) a(j,j)=tiny
          dum=1.d0/a(j,j)
          do 18 i=j+1,n
            a(i,j)=a(i,j)*dum
  18      continue
        endif
  19  continue
      if(a(n,n).eq.0.d0) a(n,n)=tiny
      return
      end
c ------------------------------------------------------------------
      subroutine lubksb(a,n,np,indx,b)
      implicit real*8(a-h,o-z)
      integer n,np,indx(np)
      real*8 a(np,np),b(n)
      ii=0
      do 12 i=1,n
        ll=indx(i)
        sum=b(ll)
        b(ll)=b(i)
        if(ii.ne.0) then
          do 11 j=ii,i-1
            sum=sum-a(i,j)*b(j)
  11      continue
        else if(sum.ne.0.d0) then
          ii=i
        endif
        b(i)=sum
  12  continue
      do 14 i=n,1,-1
        sum=b(i)
        if(i.lt.n) then
          do 13 j=i+1,n
            sum=sum-a(i,j)*b(j)
  13      continue
        endif
        b(i)=sum/a(i,i)
  14  continue
      return
      end

ext: f

max_stars_repo_name: noboruatkek/SAD

max_issues_repo_name: noboruatkek/SAD

max_forks_repo_name: noboruatkek/SAD


text:
C
C	--- DIAMOND U23 SORTING PROGRAM --- v1.91 --- B. Diviacco --- 17/8/2005
C
C	Developed and tested under Linux using the gcc compiler with f77 extension
C
C	to compile with:	g77 -fno-automatic -O u23_sort.for -o u23_sort.exe
C
C	Arrays dimensioned for:
C
C	25 X-points (+/-30 mm, step 2.5 mm)
C	Max 999 blocks of each type in the measurement data files
C	(check also that max block number (label) doensn't exceed 999)
C
C	1/8/2005:	Introduced Diamond block coordinate system (mx, mz, ms)
C			for Helmholtz Coil measurement data
C
C	2/8/2005:	Modified MOVE2 and MOVE4 to deal with D-C-B-A termination
C			Changed IQMAX from 4 to 2 in sub PMU
C			Purged a few reduntant COMMON blocks and REAL statements
C
C	9/8/2005:	Corrected block type naming and format of "setmag.out"
C
C	11/8/2005:	Added calculation of average field and rms field error
C
C	17/8/2005:	Increased array dimensions to deal with up to 999 blocks/type
C			and up to 9999 field points
C
c	COMMON/PMUDAT/X1(424,2),X2(424,2),Y1(424,2),Y2(424,2),
c     $	Z1(424,2),Z2(424,2),MX(424,2),MY(424,2),MZ(424,2)
C
	COMMON/OUT/FIRSTIX(-12:12),FIRSTIY(-12:12),
     $	SECNDIX(-12:12),SECNDIY(-12:12),BUNDX(0:9999),
     $	BUNDY(0:9999),PHQUAD(2),PHERR,FIRSTIXMAX,
     $	FIRSTIYMAX,SECNDIXMAX,SECNDIYMAX,STRX,STRY
	COMMON/PAR/PERIOD,GAP,HEIGHT,WIDTH,SEP,NP,BS1,NB
	COMMON/CONF/ITYPE(424,2),IDIR(424,2),IOR(424,2),IBLOCK(424,2),
     $	IUSE(999,4),JBLOCK(999,4),NAEXT,NBEXT,NCEXT,NDEXT
	COMMON/ZPAR/ZMIN,ZMAX,DELZ,ANEX,NSTP,NZSTEP
	COMMON/COEFF/C1,C2,C3,C4,C5,C6,C7,C8,C9,C10,C11,C12
	CHARACTER*20 FILNAM1, FILNAM2,FILNAM3
	CHARACTER*1 ORLAB
	CHARACTER*2 DIRLAB
	CHARACTER*5 TYPLAB
	LOGICAL ANS
C
C	TO EDIT: undulator geometrical parameters
C
C	BS1 is the size of the end blocks (reduced size H block).
C
	NP=90
	GAP=5.0
	PERIOD=22.0
	HEIGHT=16
	WIDTH=41
	BS1=3.7
C
C	NB (the number of blocks per array) depends on the end scheme, as defined
C	in the SETPMU subroutine.
C
	NB=4*NP+1
	WRITE(6,*)'NB=',NB
C
C	Number of steps per period and number of extra periods
C	on each side for trajectory and phase error calculation.
C
	NSTP=20
	ANEX=6
C
	open(unit=5,file='u22_I13i_sort.inp',status='old')
	READ (5,*)IAN
	IF (IAN.NE.1.AND.IAN.NE.2) STOP
	READ (5,*)NOVER, NLIMIT, NSTEP, TFACTR
	READ (5,*)C1,C2
	READ (5,*)C3,C4
	READ (5,*)C5,C6
	READ (5,*)C7,C8
	READ (5,*)C9,C10
	READ (5,*)C11,C12
	close(5)
C
C	Output file names
C
	FILNAM1='bxy.dat'
	FILNAM2='int1.dat'
	FILNAM3='int2.dat'
C
	ZMIN=-(NP+2*ANEX)*PERIOD/2
	ZMAX=-ZMIN
	NZSTEP=NSTP*(NP+2*ANEX)
	WRITE(6,*)'NZSTEP =',NZSTEP
	DELZ=(ZMAX-ZMIN)/NZSTEP
C
	CALL READ_DATA
C
	CALL READ_CONFIG
C
	CALL SETPMU
C
	CALL SETUP
C
	CALL COST(CST1)
C
	BX0=BUNDX(NZSTEP/2)
	BY0=BUNDY(NZSTEP/2)
C
C       Compute rms field error
C
	NPOLMAX=2*NP+1
	NSKIP=8
	J1=NSKIP+1
	J2=NPOLMAX-NSKIP
	NDATA=J2-J1+1
C
C	Compute Average Peak Field and rms field error
C
	BYAV=0.0
	DO J=J1,J2
	IZP=NZSTEP/2-NP*NSTP/2+(J-1)*NSTP/2
	BYAV=BYAV+ABS(BUNDY(IZP))
	END DO
	BYAV=BYAV/NDATA
C
	S=0.0
	SS=0.0
	DO J=J1,J2
	IZP=NZSTEP/2-NP*NSTP/2+(J-1)*NSTP/2
	ERRY1=ABS(BUNDY(IZP))-BYAV
	S=S+ERRY1
	SS=SS+ERRY1*ERRY1
	END DO
c	SIGMAB=SQRT(SS/NDATA-S/NDATA*S/NDATA)/BYAV*100
	SIGMAB=SQRT(SS/(NDATA-1))/BYAV*100
C
	OPEN (UNIT=66,FILE='u22_I13i_sort.out',STATUS='NEW')
	WRITE(6,*)'STARTING CONFIGURATION COST =',CST1
	WRITE(6,*)
	WRITE(6,*)'Bxo     (G)     =',BX0*1E4
	WRITE(6,*)'Byo     (G)     =',BY0*1E4
	WRITE(6,*)'<Byo>   (G)     =',BYAV*1E4
	WRITE(6,*)'Sig_B   (%)     =',SIGMAB
	WRITE(6,*)'IXmax   (G m)   =',FIRSTIXMAX*1E-3
	WRITE(6,*)'IYmax   (G m)   =',FIRSTIYMAX*1E-3
	WRITE(6,*)'IX0     (G m)   =',FIRSTIX(0)*1E-3
	WRITE(6,*)'IY0     (G m)   =',FIRSTIY(0)*1E-3
	WRITE(6,*)'I2Xmax  (G m^2) =',SECNDIXMAX*1E-6
	WRITE(6,*)'I2Ymax  (G m^2) =',SECNDIYMAX*1E-6
	WRITE(6,*)'I2X0    (G m^2) =',SECNDIX(0)*1E-6
	WRITE(6,*)'I2Y0    (G m^2) =',SECNDIY(0)*1E-6
	WRITE(6,*)'RMS Phase Error (deg) (ALL)=',PHERR
	WRITE(6,*)'RMS Phase Error (deg) (TOP)=',PHQUAD(1)
	WRITE(6,*)'RMS Phase Error (deg) (BOTTOM)=',PHQUAD(2)
	WRITE(6,*)'X-traj straightness (micron) =',STRX*1E3
	WRITE(6,*)'Y-traj straightness (micron) =',STRY*1E3
	WRITE(6,*)
C
	WRITE(66,*)'STARTING CONFIGURATION COST =',CST1
	WRITE(66,*)
	WRITE(66,*)'Bxo     (G)     =',BX0*1E4
	WRITE(66,*)'Byo     (G)     =',BY0*1E4
	WRITE(66,*)'<Byo>   (G)     =',BYAV*1E4
	WRITE(66,*)'Sig_B   (%)     =',SIGMAB
	WRITE(66,*)'IXmax   (G m)   =',FIRSTIXMAX*1E-3
	WRITE(66,*)'IYmax   (G m)   =',FIRSTIYMAX*1E-3
	WRITE(66,*)'IX0     (G m)   =',FIRSTIX(0)*1E-3
	WRITE(66,*)'IY0     (G m)   =',FIRSTIY(0)*1E-3
	WRITE(66,*)'I2Xmax  (G m^2) =',SECNDIXMAX*1E-6
	WRITE(66,*)'I2Ymax  (G m^2) =',SECNDIYMAX*1E-6
	WRITE(66,*)'I2X0    (G m^2) =',SECNDIX(0)*1E-6
	WRITE(66,*)'I2Y0    (G m^2) =',SECNDIY(0)*1E-6
	WRITE(66,*)'RMS Phase Error (deg) (ALL)=',PHERR
	WRITE(66,*)'RMS Phase Error (deg) (TOP)=',PHQUAD(1)
	WRITE(66,*)'RMS Phase Error (deg) (BOTTOM)=',PHQUAD(2)
	WRITE(66,*)'X-traj straightness (micron) =',STRX*1E3
	WRITE(66,*)'Y-traj straightness (micron) =',STRY*1E3
	WRITE(66,*)
C
	IF (IAN.EQ.2) GOTO 888
C
	TSTART=1.0E6
C
	WRITE(6,*)
	WRITE(6,*)'*** STARTING ANNEALING PROCESS ***'
	WRITE(6,*)
C
	OPEN(UNIT=3,FILE='history.dat',STATUS='NEW')
C
C	Annealing schedule
C
	T=TSTART
C	R=RAND(1234)
	I = 86456
	I=time()
	CALL SRAND(I)
	I=0
C
C	Try up to NSTEP temperature steps...
C
	NTOT=1
	DO 6 ISTEP=1,NSTEP
	NSUCC=0
	NSUCC1=0
	NSUCC2=0
	NSUCC3=0
	NSUCC4=0
	NSUCC5=0
	NSUCC6=0
	AVCOST=0.0
	NTRY=0
	DEMEAN=0.0
C
C	... and up to NOVER reconfigurations.
C
	DO 4 NREC=1,NOVER
C
C 	Propose a Random Move
C
777	IMOVE=INT(RAND(0)*6.0)+1
	IF (IMOVE.EQ.1) CALL MOVE1(1)
	IF (IMOVE.EQ.2) CALL MOVE2(1)
	IF (IMOVE.EQ.3) CALL MOVE3(1)
	IF (IMOVE.EQ.4) CALL MOVE4(1)
	IF (IMOVE.EQ.5) CALL MOVE5(1)
	IF (IMOVE.EQ.6) CALL MOVE6(1)
C
	CALL COST(CST2)

	NTOT=NTOT+1
C
C	Consult the oracle
C
	DE=CST2-CST1
	CALL METROP(DE,T,ANS)
C
	IF (ANS) THEN
C
C	If move is accepted then...
C
	IF (IMOVE.EQ.1) NSUCC1=NSUCC1+1
	IF (IMOVE.EQ.2) NSUCC2=NSUCC2+1
	IF (IMOVE.EQ.3) NSUCC3=NSUCC3+1
	IF (IMOVE.EQ.4) NSUCC4=NSUCC4+1
	IF (IMOVE.EQ.5) NSUCC5=NSUCC5+1
	IF (IMOVE.EQ.6) NSUCC6=NSUCC6+1
	NSUCC=NSUCC1+NSUCC2+NSUCC3+NSUCC4+NSUCC5+NSUCC6
	AVCOST=AVCOST+CST2
	NTRY=NTRY+1
	DEMEAN=DEMEAN+ABS(DE)
	CST1=CST2
	ELSE
C
C	...else go back to the previous configuration.
C
	IF (IMOVE.EQ.1) CALL MOVE1(0)
	IF (IMOVE.EQ.2) CALL MOVE2(0)
	IF (IMOVE.EQ.3) CALL MOVE3(0)
	IF (IMOVE.EQ.4) CALL MOVE4(0)
	IF (IMOVE.EQ.5) CALL MOVE5(0)
	IF (IMOVE.EQ.6) CALL MOVE6(0)
	ENDIF
C
	IF (NSUCC.GE.NLIMIT) GOTO 5
4	CONTINUE
5	IF(NTRY.NE.0)DEMEAN=DEMEAN/NTRY
	IF(NTRY.NE.0)AVCOST=AVCOST/NTRY
	WRITE(6,110)ISTEP,'T = ',T,'       FINAL COST = ',CST1
	WRITE(6,112)'      AVERAGE COST =',AVCOST
	WRITE(6,103)DEMEAN
	WRITE(6,111)'Successful M1 moves: ',NSUCC1
	WRITE(6,111)'Successful M2 moves: ',NSUCC2
	WRITE(6,111)'Successful M3 moves: ',NSUCC3
	WRITE(6,111)'Successful M4 moves: ',NSUCC4
	WRITE(6,111)'Successful M5 moves: ',NSUCC5
	WRITE(6,111)'Successful M6 moves: ',NSUCC6
	WRITE(6,111)'---------------------',NSUCC
	IF(ISTEP.NE.1)WRITE(3,*)T,AVCOST
C
C	Check if the final goal is reached
C
	IF (NSUCC.EQ.0) GOTO 7
C
	T=T*TFACTR
	IF (ISTEP.EQ.1) T=DEMEAN*10
6	CONTINUE
C
7	WRITE(6,*)
        WRITE(6,*)'***************************************************'
	WRITE(6,*)' Number of configurations examined =',NTOT
	WRITE(6,*)
	WRITE(6,*)'                  *** SYSTEM FROZEN ***'
	WRITE(6,*)
        WRITE(6,*)'***************************************************'
	WRITE(6,*)
C
	CLOSE(3)
C
	CALL COST(CST3)
C
C	WRITE CONFIGURATION
C
	OPEN(UNIT=2,FILE='setmag.out',STATUS='NEW')
	DO IQ=1,2
	DO IPOS=1,NB
	IITYPE=ITYPE(IPOS,IQ)
	IIDIR=IDIR(IPOS,IQ)
	IIOR=IOR(IPOS,IQ)
	IIBLOCK=IBLOCK(IPOS,IQ)
C
	IF (IIOR.EQ.1)  ORLAB='N'
	IF (IIOR.EQ.-1) ORLAB='R'
C
	IF (IITYPE.EQ.1) TYPLAB='Vfull'
	IF (IITYPE.EQ.2) TYPLAB='Hfull'
	IF (IITYPE.EQ.3) TYPLAB='Vhalf'
	IF (IITYPE.EQ.4) TYPLAB='Hhalf'
C
	IF ((IITYPE.EQ.1.OR.IITYPE.EQ.3).AND.IIDIR.EQ.1) DIRLAB='UP'
	IF ((IITYPE.EQ.1.OR.IITYPE.EQ.3).AND.IIDIR.EQ.-1) DIRLAB='DW'
	IF ((IITYPE.EQ.2.OR.IITYPE.EQ.4).AND.IIDIR.EQ.1) DIRLAB='RG'
	IF ((IITYPE.EQ.2.OR.IITYPE.EQ.4).AND.IIDIR.EQ.-1) DIRLAB='LF'
C
	WRITE(2,1111)IQ,IPOS,IITYPE,IIDIR,IIOR,IIBLOCK,TYPLAB,DIRLAB,ORLAB
c	WRITE(2,1112)IQ,IPOS,IITYPE,IIDIR,IIOR,IIBLOCK,TYPLAB,DIRLAB,ORLAB,
c     $	MX(IPOS,IQ),MY(IPOS,IQ),MZ(IPOS,IQ)
	END DO
	WRITE(2,*)
	END DO
	CLOSE(2)
1111	FORMAT(6(2X,I3),6X,A5,2X,A2,2X,A1)
c1112	FORMAT(6(2X,I3),6X,A5,2X,A2,2X,A1,3(2X,F6.3))
C
C	WRITE FIELD IN 'FILNAM1' AND FIELD INTEGRALS IN 'FILNAM2 AND FILNAM3'
C
888	WRITE(6,*)'WRITING FIELD AND FIELD INTEGRALS'
	OPEN(UNIT=1,FILE=FILNAM1,STATUS='NEW')
	DO IZ=0,NZSTEP
	PZ=ZMIN+(ZMAX-ZMIN)*IZ/NZSTEP
	BX=BUNDX(IZ)
	BY=BUNDY(IZ)
	WRITE(1,*)BX*1E4,BY*1E4,PZ
	END DO
	CLOSE(1)
C
	OPEN(UNIT=2,FILE=FILNAM2,STATUS='NEW')
	OPEN(UNIT=3,FILE=FILNAM3,STATUS='NEW')
	DO IX=-12,12,1
	XX=IX*2.5
	FX=FIRSTIX(IX)*1E-3
	FY=FIRSTIY(IX)*1E-3
	SX=SECNDIX(IX)*1E-6
	SY=SECNDIY(IX)*1E-6
	WRITE(2,220)XX,FX,FY
	WRITE(3,220)XX,SX,SY
	END DO
	CLOSE(2)
	CLOSE(3)
C
	IF (IAN.EQ.2) GOTO 999
C
	BX0=BUNDX(NZSTEP/2)
	BY0=BUNDY(NZSTEP/2)
C
C	Compute Average Peak Field and rms field error
C
	BYAV=0.0
	DO J=J1,J2
	IZP=NZSTEP/2-NP*NSTP/2+(J-1)*NSTP/2
	BYAV=BYAV+ABS(BUNDY(IZP))
	END DO
	BYAV=BYAV/NDATA
C
	S=0.0
	SS=0.0
	DO J=J1,J2
	IZP=NZSTEP/2-NP*NSTP/2+(J-1)*NSTP/2
	ERRY1=ABS(BUNDY(IZP))-BYAV
	S=S+ERRY1
	SS=SS+ERRY1*ERRY1
	END DO
c	SIGMAB=SQRT(SS/NDATA-S/NDATA*S/NDATA)/BYAV*100
	SIGMAB=SQRT(SS/(NDATA-1))/BYAV*100
C
	WRITE(6,*)
	WRITE(6,*)'FINAL CONFIGURATION COST =',CST3
	WRITE(6,*)
	WRITE(6,*)'Bxo     (G)     =',BX0*1E4
	WRITE(6,*)'Byo     (G)     =',BY0*1E4
	WRITE(6,*)'<Byo>   (G)     =',BYAV*1E4
	WRITE(6,*)'Sig_B   (%)     =',SIGMAB
	WRITE(6,*)'IXmax   (G m)  =',FIRSTIXMAX*1E-3,'( C1=',C1,')'
	WRITE(6,*)'IYmax   (G m)  =',FIRSTIYMAX*1E-3,'( C2=',C2,')'
	WRITE(6,*)'IX0     (G m)  =',FIRSTIX(0)*1E-3,'( C3=',C3,')'
	WRITE(6,*)'IY0     (G m)  =',FIRSTIY(0)*1E-3,'( C4=',C4,')'
	WRITE(6,*)'I2Xmax  (G m^2) =',SECNDIXMAX*1E-6,'( C5=',C5,')'
	WRITE(6,*)'I2Ymax  (G m^2) =',SECNDIYMAX*1E-6,'( C6=',C6,')'
	WRITE(6,*)'I2X0    (G m^2) =',SECNDIX(0)*1E-6,'( C7=',C7,')'
	WRITE(6,*)'I2Y0    (G m^2) =',SECNDIY(0)*1E-6,'( C8=',C8,')'
	WRITE(6,*)'RMS Phase Error (deg) (ALL)=',PHERR,'( C9=',C9,')'
	WRITE(6,*)'RMS Phase Error (deg) (TOP)=',PHQUAD(1),'( C10=',C10,')'
	WRITE(6,*)'RMS Phase Error (deg) (BOTTOM)=',PHQUAD(2)
	WRITE(6,*)'X-traj straightness (micron) =',STRX*1E3,'( C11=',C11,')'
	WRITE(6,*)'Y-traj straightness (micron) =',STRY*1E3,'( C12=',C12,')'
	WRITE(6,*)
C
	WRITE(66,*)
	WRITE(66,*)'FINAL CONFIGURATION COST =',CST3
	WRITE(66,*)
	WRITE(66,*)'Bxo     (G)     =',BX0*1E4
	WRITE(66,*)'Byo     (G)     =',BY0*1E4
	WRITE(66,*)'<Byo>   (G)     =',BYAV*1E4
	WRITE(66,*)'Sig_B   (%)     =',SIGMAB
	WRITE(66,*)'IXmax   (G m)  =',FIRSTIXMAX*1E-3,'( C1=',C1,')'
	WRITE(66,*)'IYmax   (G m)  =',FIRSTIYMAX*1E-3,'( C2=',C2,')'
	WRITE(66,*)'IX0     (G m)  =',FIRSTIX(0)*1E-3,'( C3=',C3,')'
	WRITE(66,*)'IY0     (G m)  =',FIRSTIY(0)*1E-3,'( C4=',C4,')'
	WRITE(66,*)'I2Xmax  (G m^2) =',SECNDIXMAX*1E-6,'( C5=',C5,')'
	WRITE(66,*)'I2Ymax  (G m^2) =',SECNDIYMAX*1E-6,'( C6=',C6,')'
	WRITE(66,*)'I2X0    (G m^2) =',SECNDIX(0)*1E-6,'( C7=',C7,')'
	WRITE(66,*)'I2Y0    (G m^2) =',SECNDIY(0)*1E-6,'( C8=',C8,')'
	WRITE(66,*)'RMS Phase Error (deg) (ALL)=',PHERR,'( C9=',C9,')'
	WRITE(66,*)'RMS Phase Error (deg) (TOP)=',PHQUAD(1),'( C10=',C10,')'
	WRITE(66,*)'RMS Phase Error (deg) (BOTTOM)=',PHQUAD(2)
	WRITE(66,*)'X-traj straightness (micron) =',STRX*1E3,'( C11=',C11,')'
	WRITE(66,*)'Y-traj straightness (micron) =',STRY*1E3,'( C12=',C12,')'
	WRITE(66,*)
C
999	CLOSE(66)
	STOP
C
103	FORMAT(34X,'Typical Delta_E =',F12.6)
110     FORMAT(1X,I3,2X,A4,F12.4,10X,A20,F15.6)
111     FORMAT(1X,A21,I6,A6,I6)
112     FORMAT(32X,A20,F15.6)
220	FORMAT(3(2X,F12.3))
	END
C
C	--------------------
	SUBROUTINE READ_DATA
C	--------------------
	REAL MMX,MMY,MMZ
	COMMON/CONF/ITYPE(424,2),IDIR(424,2),IOR(424,2),IBLOCK(424,2),
     $	IUSE(999,4),JBLOCK(999,4),NAEXT,NBEXT,NCEXT,NDEXT
	COMMON/ARCHIVE/MMX(999,4),MMY(999,4),MMZ(999,4)
	COMMON/PAR/PERIOD,GAP,HEIGHT,WIDTH,SEP,NP,BS1,NB

	COMMON/BLOCKNUM/NAMIN,NBMIN,NCMIN,NDMIN,NAMAX,NBMAX,NCMAX,NDMAX,
     $	NBA,NBB,NBC,NBD
C
	WRITE(6,*)'READING MEASUREMENT DATA FILES'
C
C	TO EDIT: Total number of blocks for each type in the measurement database
C
C	A=vertically magnetized, full size (V-full)
C	B=horizontally magnetized, full size (H-full)
C	C=vertically magnetized, half size (V-half)
C	D=horizontally magnetized, reduced size (H-half)
C
	NBA=424
	NBB=424
	NBC=12
	NBD=12
C
C	TO EDIT: Minimum and Maximum block # for each block type
C
	NAMIN=1
	NBMIN=3
	NCMIN=1
	NDMIN=2
C
	NAMAX=466
	NBMAX=677
	NCMAX=12
	NDMAX=21
C
	DO I=NAMIN,NAMAX
	IUSE(I,1)=0
	END DO
C
	DO I=NBMIN,NBMAX
	IUSE(I,2)=0
	END DO
C
	DO I=NCMIN,NCMAX
	IUSE(I,3)=0
	END DO
C
	DO I=NDMIN,NDMAX
	IUSE(I,4)=0
	END DO
C
C	TO EDIT: magnetic data filenames and/or fileformat
C
	OPEN(UNIT=1,FILE='I13iV.sim',STATUS='OLD')
	DO I=1,NBA
	READ(1,*)DUM,AX,AZ,AS
	N=INT(DUM)
	MMX(N,1)=AX
	MMY(N,1)=AZ
	MMZ(N,1)=AS
	IUSE(N,1)=2
	END DO
	CLOSE(1)
C
	OPEN(UNIT=2,FILE='I13iH.sim',STATUS='OLD')
	DO I=1,NBB
	READ(2,*)DUM,AX,AZ,AS
	N=INT(DUM)
	MMX(N,2)=AX
	MMY(N,2)=AZ
	MMZ(N,2)=AS
	IUSE(N,2)=2
	END DO
	CLOSE(2)
C
	OPEN(UNIT=3,FILE='I13iVE.sim',STATUS='OLD')
	DO I=1,NBC
	READ(3,*)DUM,AX,AZ,AS
	N=INT(DUM)
	MMX(N,3)=AX
	MMY(N,3)=AZ
	MMZ(N,3)=AS
	IUSE(N,3)=2
	END DO
	CLOSE(3)
C
	OPEN(UNIT=4,FILE='I13iHEA.sim',STATUS='OLD')
	DO I=1,NBD
	READ(4,*)DUM,AX,AZ,AS
	N=INT(DUM)
	MMX(N,4)=AX
	MMY(N,4)=AZ
	MMZ(N,4)=AS
	IUSE(N,4)=2
	END DO
	CLOSE(4)
C
C	TO EDIT: Mark blocks which should not be used, even if
C	present in the measurement database (for example damaged
C	or incorrectly measured blocks).
C
C	First argument of IUSE is the block number (label)
C	Second index is the block type (A=1, B=2, C=3, D=4)
C
c	IUSE(3,1)=0
c	IUSE(97,2)=0
C
	RETURN
	END
C
C	----------------------
	SUBROUTINE READ_CONFIG
C	----------------------
	COMMON/PAR/PERIOD,GAP,HEIGHT,WIDTH,SEP,NP,BS1,NB
	COMMON/CONF/ITYPE(424,2),IDIR(424,2),IOR(424,2),IBLOCK(424,2),
     $	IUSE(999,4),JBLOCK(999,4),NAEXT,NBEXT,NCEXT,NDEXT
	COMMON/BLOCKNUM/NAMIN,NBMIN,NCMIN,NDMIN,NAMAX,NBMAX,NCMAX,NDMAX,
     $	NBA,NBB,NBC,NBD
	LOGICAL L1,L2,L3,L4
C
	WRITE(6,*)'READING setmag.inp'
	OPEN(UNIT=1,FILE='setmag.inp',STATUS='OLD')
	DO IQ=1,2
	DO IPOS=1,NB
	READ (1,111)IIQ,IIPOS,IITYPE,IIDIR,IIOR,IIBLOCK
	IF (IIQ.NE.IQ.OR.IIPOS.NE.IPOS) GOTO 444
	IF (IUSE(IIBLOCK,IITYPE).EQ.0) GOTO 666
	ITYPE(IPOS,IQ)=IITYPE
	IDIR(IPOS,IQ)=IIDIR
	IOR(IPOS,IQ)=IIOR
	IBLOCK(IPOS,IQ)=IIBLOCK
C
C	Mark blocks used in the initial configuration
C
	IUSE(IIBLOCK,IITYPE)=1
c
	END DO
	READ (1,*)
	END DO
	CLOSE(1)
111	FORMAT(6(2X,I3))
C
C	CHECK BLOCK NUMBERING
C
	WRITE(6,*) 'CHECKING UNICITY OF BLOCK NUMBERING'
	WRITE(6,*)
	DO IPOS1=1,NB
	DO IQ1=1,2
	IBLOCK1=IBLOCK(IPOS1,IQ1)
	ITYPE1=ITYPE(IPOS1,IQ1)
	DO IPOS2=1,NB
	DO IQ2=1,2
	IBLOCK2=IBLOCK(IPOS2,IQ2)
	ITYPE2=ITYPE(IPOS2,IQ2)
	L1=(IBLOCK2.EQ.IBLOCK1)
	L2=(ITYPE2.EQ.ITYPE1)
	L3=(IQ2.NE.IQ1)
	L4=(IPOS2.NE.IPOS1)
	IF (L1.AND.L2.AND.(L3.OR.L4)) GOTO 555
	END DO
	END DO
	END DO
	END DO
c
	nanot=0
	nbnot=0
	ncnot=0
	ndnot=0
	naconf=0
	nbconf=0
	ncconf=0
	ndconf=0
	naext=0
	nbext=0
	ncext=0
	ndext=0
	do i=namin,namax
	if (iuse(i,1).eq.0) nanot=nanot+1
	if (iuse(i,1).eq.1) naconf=naconf+1
	if (iuse(i,1).eq.2) then
	naext=naext+1
	jblock(naext,1)=i
	endif
	end do
	do i=nbmin,nbmax
	if (iuse(i,2).eq.0) nbnot=nbnot+1
	if (iuse(i,2).eq.1) nbconf=nbconf+1
	if (iuse(i,2).eq.2) then
	nbext=nbext+1
	jblock(nbext,2)=i
	endif
	end do
	do i=ncmin,ncmax
	if (iuse(i,3).eq.0) ncnot=ncnot+1
	if (iuse(i,3).eq.1) ncconf=ncconf+1
	if (iuse(i,3).eq.2) then
	ncext=ncext+1
	jblock(ncext,3)=i
	endif
	end do
	do i=ndmin,ndmax
	if (iuse(i,4).eq.0) ndnot=ndnot+1
	if (iuse(i,4).eq.1) ndconf=ndconf+1
	if (iuse(i,4).eq.2) then
	ndext=ndext+1
	jblock(ndext,4)=i
	endif
	end do
	write(6,*)'A blocks not usable or missing',nanot
	write(6,*)'B blocks not usable or missing',nbnot
	write(6,*)'C blocks not usable or missing',ncnot
	write(6,*)'D blocks not usable or missing',ndnot
	write(6,*)
	write(6,*)'A blocks in the configuration',naconf
	write(6,*)'B blocks in the configuration',nbconf
	write(6,*)'C blocks in the configuration',ncconf
	write(6,*)'D blocks in the configuration',ndconf
	write(6,*)
	write(6,*)'A blocks available for swapping',naext
	write(6,*)'B blocks available for swapping',nbext
	write(6,*)'C blocks available for swapping',ncext
	write(6,*)'D blocks available for swapping',ndext
	write(6,*)
C
	RETURN
444	WRITE(6,*)'ERROR WHILE READING setmag.inp'
	STOP
555	WRITE(6,*)'*** BLOCK NUMBERING ERROR ***'
	WRITE(6,*)' *  block at array ',iq1,'    ,position ',ipos1,
     $	'   type ',itype1
	WRITE(6,*)' *  has same label (',itype1,') as'
	WRITE(6,*)' *  block at array ',iq2,'    ,position ',ipos2,
     $	'   type ',itype2
	STOP
666	WRITE(6,*)'*** BLOCK LABEL ERROR ***'
	WRITE(6,*)' *  block ',iiblock,' itype ',iitype,
     $	' should not be used !'
	STOP
	END
C
C	-----------------
	SUBROUTINE SETPMU
C	-----------------
	REAL MX,MY,MZ
	COMMON/PMUDAT/X1(424,2),X2(424,2),Y1(424,2),Y2(424,2),
     $	Z1(424,2),Z2(424,2),MX(424,2),MY(424,2),MZ(424,2)
	COMMON/PAR/PERIOD,GAP,HEIGHT,WIDTH,SEP,NP,BS1,NB
	COMMON/CONF/ITYPE(424,2),IDIR(424,2),IOR(424,2),IBLOCK(424,2),
     $	IUSE(999,4),JBLOCK(999,4),NAEXT,NBEXT,NCEXT,NDEXT
C
C	SET UP BLOCK COORDINATES IN COMMON/PMUDAT/
C
C	WRITE(6,*)'SETTING UP BLOCK COORDINATES'
C
	W=(PERIOD/4)
C
	DO IQ=1,2
	ZMIN=-((NB*W)/2)
	XMIN=-WIDTH/2
	XMAX=WIDTH/2
	IF (IQ.EQ.1) THEN
	YMIN=GAP/2.0
	YMAX=YMIN+HEIGHT
	ELSE
	YMAX=-GAP/2.0
	YMIN=YMAX-HEIGHT
	ENDIF
C
	DO I=1,NB
	IF(I.EQ.1) THEN
	Z1(I,IQ)=ZMIN
	Z2(I,IQ)=ZMIN+W
	ELSE
	Z1(I,IQ)=Z2(I-1,IQ)
	Z2(I,IQ)=Z1(I,IQ)+W
	ENDIF
C
	Y1(I,IQ)=YMIN
	Y2(I,IQ)=YMAX
	X1(I,IQ)=XMIN
	X2(I,IQ)=XMAX
C
	END DO
C
	I1=1
	I2=2
	I3=NB-1
	I4=NB
C
C	Adjust size and position of the end-blocks
C
	Z1(I2,IQ)=Z2(I2,IQ)-W/2
	Z2(I1,IQ)=Z1(I2,IQ)
	Z1(I1,IQ)=Z2(I1,IQ)-BS1
C
	Z2(I3,IQ)=Z1(I3,IQ)+W/2
	Z1(I4,IQ)=Z2(I3,IQ)
	Z2(I4,IQ)=Z1(I4,IQ)+BS1
C
	END DO
C
	RETURN
	END
C
C	----------------
	SUBROUTINE SETUP
C	----------------
	REAL MX,MY,MZ
	COMMON/PMUDAT/X1(424,2),X2(424,2),Y1(424,2),Y2(424,2),
     $	Z1(424,2),Z2(424,2),MX(424,2),MY(424,2),MZ(424,2)
	COMMON/PAR/PERIOD,GAP,HEIGHT,WIDTH,SEP,NP,BS1,NB
	COMMON/CONF/ITYPE(424,2),IDIR(424,2),IOR(424,2),IBLOCK(424,2),
     $	IUSE(999,4),JBLOCK(999,4),NAEXT,NBEXT,NCEXT,NDEXT
	COMMON/ZPAR/ZMIN,ZMAX,DELZ,ANEX,NSTP,NZSTEP
	COMMON/FLD/BBX(424,2,3,0:9999),BBY(424,2,3,0:9999),
     $	FINTX(424,2,3,-12:12),FINTY(424,2,3,-12:12),
     $	SINTX(424,2,3,-12:12),SINTY(424,2,3,-12:12)
	DIMENSION C(3)
C
	WRITE(6,*)
	WRITE(6,*)'COMPUTING BLOCK SIGNATURES'
C
	C(2)=0.0
C
	DO IQ=1,2
	WRITE(6,*)'ARRAY =',IQ
	DO I=1,NB
C
	DO IX=-12,12,1
	C(1)=IX*2.5
C
	FINTX(I,IQ,1,IX)=0.0
	FINTX(I,IQ,2,IX)=0.0
	FINTX(I,IQ,3,IX)=0.0
	FINTY(I,IQ,1,IX)=0.0
	FINTY(I,IQ,2,IX)=0.0
	FINTY(I,IQ,3,IX)=0.0
C
	SINTX(I,IQ,1,IX)=0.0
	SINTX(I,IQ,2,IX)=0.0
	SINTX(I,IQ,3,IX)=0.0
	SINTY(I,IQ,1,IX)=0.0
	SINTY(I,IQ,2,IX)=0.0
	SINTY(I,IQ,3,IX)=0.0
C
	DO IZ=0,NZSTEP
	C(3)=ZMIN+(ZMAX-ZMIN)*IZ/NZSTEP
	MX(I,IQ)=1.0
	MY(I,IQ)=0.0
	MZ(I,IQ)=0.0
	CALL PMU(C,BX,1,I,IQ)
	CALL PMU(C,BY,2,I,IQ)
	IF (IX.EQ.0) THEN
	BBX(I,IQ,1,IZ)=BX
	BBY(I,IQ,1,IZ)=BY
	ENDIF
	FINTX(I,IQ,1,IX)=FINTX(I,IQ,1,IX)+BX*1E4*DELZ
	FINTY(I,IQ,1,IX)=FINTY(I,IQ,1,IX)+BY*1E4*DELZ
	SINTX(I,IQ,1,IX)=SINTX(I,IQ,1,IX)-BX*1E4*DELZ*C(3)
	SINTY(I,IQ,1,IX)=SINTY(I,IQ,1,IX)-BY*1E4*DELZ*C(3)
C
	MX(I,IQ)=0.0
	MY(I,IQ)=1.0
	MZ(I,IQ)=0.0
	CALL PMU(C,BX,1,I,IQ)
	CALL PMU(C,BY,2,I,IQ)
	IF (IX.EQ.0) THEN
	BBX(I,IQ,2,IZ)=BX
	BBY(I,IQ,2,IZ)=BY
	ENDIF
	FINTX(I,IQ,2,IX)=FINTX(I,IQ,2,IX)+BX*1E4*DELZ
	FINTY(I,IQ,2,IX)=FINTY(I,IQ,2,IX)+BY*1E4*DELZ
	SINTX(I,IQ,2,IX)=SINTX(I,IQ,2,IX)-BX*1E4*DELZ*C(3)
	SINTY(I,IQ,2,IX)=SINTY(I,IQ,2,IX)-BY*1E4*DELZ*C(3)
C
	MX(I,IQ)=0.0
	MY(I,IQ)=0.0
	MZ(I,IQ)=1.0
	CALL PMU(C,BX,1,I,IQ)
	CALL PMU(C,BY,2,I,IQ)
	IF (IX.EQ.0) THEN
	BBX(I,IQ,3,IZ)=BX
	BBY(I,IQ,3,IZ)=BY
	ENDIF
	FINTX(I,IQ,3,IX)=FINTX(I,IQ,3,IX)+BX*1E4*DELZ
	FINTY(I,IQ,3,IX)=FINTY(I,IQ,3,IX)+BY*1E4*DELZ
	SINTX(I,IQ,3,IX)=SINTX(I,IQ,3,IX)-BX*1E4*DELZ*C(3)
	SINTY(I,IQ,3,IX)=SINTY(I,IQ,3,IX)-BY*1E4*DELZ*C(3)
C
	END DO
	END DO
	END DO
	END DO
C
	RETURN
	END
C
C	-----------------------
	SUBROUTINE COST(CST)
C	-----------------------
	REAL MX,MY,MZ
	COMMON/PAR/PERIOD,GAP,HEIGHT,WIDTH,SEP,NP,BS1,NB
	COMMON/PMUDAT/X1(424,2),X2(424,2),Y1(424,2),Y2(424,2),
     $	Z1(424,2),Z2(424,2),MX(424,2),MY(424,2),MZ(424,2)
	COMMON/CONF/ITYPE(424,2),IDIR(424,2),IOR(424,2),IBLOCK(424,2),
     $	IUSE(999,4),JBLOCK(999,4),NAEXT,NBEXT,NCEXT,NDEXT
	COMMON/ZPAR/ZMIN,ZMAX,DELZ,ANEX,NSTP,NZSTEP
	COMMON/FLD/BBX(424,2,3,0:9999),BBY(424,2,3,0:9999),
     $	FINTX(424,2,3,-12:12),FINTY(424,2,3,-12:12),
     $	SINTX(424,2,3,-12:12),SINTY(424,2,3,-12:12)
	COMMON/OUT/FIRSTIX(-12:12),FIRSTIY(-12:12),
     $	SECNDIX(-12:12),SECNDIY(-12:12),BUNDX(0:9999),
     $	BUNDY(0:9999),PHQUAD(2),PHERR,FIRSTIXMAX,
     $	FIRSTIYMAX,SECNDIXMAX,SECNDIYMAX,STRX,STRY
	COMMON/TRAJ/X(0:9999),XP(0:9999),Y(0:9999),YP(0:9999)
	COMMON/COEFF/C1,C2,C3,C4,C5,C6,C7,C8,C9,C10,C11,C12
C
	CALL ASSIGN
C
C	Individual arrays
C
	DO IQ=1,2
	DO IZ=0,NZSTEP
	BX=0.0
	BY=0.0
	DO IPOS=1,NB
	BX=BX+BBX(IPOS,IQ,1,IZ)*MX(IPOS,IQ)
	BX=BX+BBX(IPOS,IQ,2,IZ)*MY(IPOS,IQ)
	BX=BX+BBX(IPOS,IQ,3,IZ)*MZ(IPOS,IQ)
	BY=BY+BBY(IPOS,IQ,1,IZ)*MX(IPOS,IQ)
	BY=BY+BBY(IPOS,IQ,2,IZ)*MY(IPOS,IQ)
	BY=BY+BBY(IPOS,IQ,3,IZ)*MZ(IPOS,IQ)
	END DO
	BUNDX(IZ)=BX
	BUNDY(IZ)=BY
	END DO
C
	CALL PHASE
	PHQUAD(IQ)=PHERR
	END DO

C	Complete Undulator
C
	DO IX=-12,12,1
C
	FIRSTIX(IX)=0.0
	FIRSTIY(IX)=0.0
	SECNDIX(IX)=0.0
	SECNDIY(IX)=0.0
C
	IF (IX.EQ.0) THEN
C
	DO IZ=0,NZSTEP
	BX=0.0
	BY=0.0
	DO IQ=1,2
	DO IPOS=1,NB
	BX=BX+BBX(IPOS,IQ,1,IZ)*MX(IPOS,IQ)
	BX=BX+BBX(IPOS,IQ,2,IZ)*MY(IPOS,IQ)
	BX=BX+BBX(IPOS,IQ,3,IZ)*MZ(IPOS,IQ)
	BY=BY+BBY(IPOS,IQ,1,IZ)*MX(IPOS,IQ)
	BY=BY+BBY(IPOS,IQ,2,IZ)*MY(IPOS,IQ)
	BY=BY+BBY(IPOS,IQ,3,IZ)*MZ(IPOS,IQ)
	END DO
	END DO
	BUNDX(IZ)=BX
	BUNDY(IZ)=BY
	END DO
C
	ENDIF
C
	DO IQ=1,2
	DO IPOS=1,NB
	FIRSTIX(IX)=FIRSTIX(IX)+FINTX(IPOS,IQ,1,IX)*MX(IPOS,IQ)
	FIRSTIX(IX)=FIRSTIX(IX)+FINTX(IPOS,IQ,2,IX)*MY(IPOS,IQ)
	FIRSTIX(IX)=FIRSTIX(IX)+FINTX(IPOS,IQ,3,IX)*MZ(IPOS,IQ)
	FIRSTIY(IX)=FIRSTIY(IX)+FINTY(IPOS,IQ,1,IX)*MX(IPOS,IQ)
	FIRSTIY(IX)=FIRSTIY(IX)+FINTY(IPOS,IQ,2,IX)*MY(IPOS,IQ)
	FIRSTIY(IX)=FIRSTIY(IX)+FINTY(IPOS,IQ,3,IX)*MZ(IPOS,IQ)
C
	SECNDIX(IX)=SECNDIX(IX)+SINTX(IPOS,IQ,1,IX)*MX(IPOS,IQ)
	SECNDIX(IX)=SECNDIX(IX)+SINTX(IPOS,IQ,2,IX)*MY(IPOS,IQ)
	SECNDIX(IX)=SECNDIX(IX)+SINTX(IPOS,IQ,3,IX)*MZ(IPOS,IQ)
	SECNDIY(IX)=SECNDIY(IX)+SINTY(IPOS,IQ,1,IX)*MX(IPOS,IQ)
	SECNDIY(IX)=SECNDIY(IX)+SINTY(IPOS,IQ,2,IX)*MY(IPOS,IQ)
	SECNDIY(IX)=SECNDIY(IX)+SINTY(IPOS,IQ,3,IX)*MZ(IPOS,IQ)
	END DO
	END DO
C
	END DO
C
	FIRSTIXMAX=0.0
	FIRSTIYMAX=0.0
	DO IX=-12,12,1
	IF (ABS(FIRSTIX(IX)).GT.FIRSTIXMAX) FIRSTIXMAX=ABS(FIRSTIX(IX))
	IF (ABS(FIRSTIY(IX)).GT.FIRSTIYMAX) FIRSTIYMAX=ABS(FIRSTIY(IX))
	END DO
	SECNDIXMAX=0.0
	SECNDIYMAX=0.0
	DO IX=-12,12,1
	IF (ABS(SECNDIX(IX)).GT.SECNDIXMAX) SECNDIXMAX=ABS(SECNDIX(IX))
	IF (ABS(SECNDIY(IX)).GT.SECNDIYMAX) SECNDIYMAX=ABS(SECNDIY(IX))
	END DO
C
	CALL PHASE
C
	CALL STRAIGHT
C
C	Cost Function definition
C
	CST=C1*FIRSTIXMAX*1E-3
	CST=CST+C2*FIRSTIYMAX*1E-3
	CST=CST+C3*ABS(FIRSTIX(0))*1E-3
	CST=CST+C4*ABS(FIRSTIY(0))*1E-3
	CST=CST+C5*SECNDIXMAX*1E-6
	CST=CST+C6*SECNDIYMAX*1E-6
	CST=CST+C7*ABS(SECNDIX(0))*1E-6
	CST=CST+C8*ABS(SECNDIY(0))*1E-6
	CST=CST+C9*PHERR
	DO IQ=1,2
	CST=CST+C10*PHQUAD(IQ)
	END DO
	CST=CST+C11*STRX*1E3+C12*STRY*1E3
C
	RETURN
	END

C
C	-----------------
	SUBROUTINE ASSIGN
C	-----------------
	REAL MX,MY,MZ,MMX,MMY,MMZ
	COMMON/PMUDAT/X1(424,2),X2(424,2),Y1(424,2),Y2(424,2),
     $	Z1(424,2),Z2(424,2),MX(424,2),MY(424,2),MZ(424,2)
	COMMON/PAR/PERIOD,GAP,HEIGHT,WIDTH,SEP,NP,BS1,NB
	COMMON/CONF/ITYPE(424,2),IDIR(424,2),IOR(424,2),IBLOCK(424,2),
     $	IUSE(999,4),JBLOCK(999,4),NAEXT,NBEXT,NCEXT,NDEXT
	COMMON/ARCHIVE/MMX(999,4),MMY(999,4),MMZ(999,4)
C
C	SET UP BLOCK MAGNETISATIONS IN COMMON/PMUDAT/
C
C	WRITE(6,*)'ASSIGNING MAGMETIZATION VECTORS'
C
	DO IQ=1,2
C
	DO I=1,NB
C
C	ASSIGN MAGNETIZATION DATA
C
	IF (ITYPE(I,IQ).EQ.1.OR.ITYPE(I,IQ).EQ.3) THEN
C
	IF (IDIR(I,IQ).EQ.1) THEN
C
	MX(I,IQ)=MMX(IBLOCK(I,IQ),ITYPE(I,IQ))
	MY(I,IQ)=MMY(IBLOCK(I,IQ),ITYPE(I,IQ))
	MZ(I,IQ)=MMZ(IBLOCK(I,IQ),ITYPE(I,IQ))
C
	ENDIF
C
	IF (IDIR(I,IQ).EQ.-1) THEN
C
	MX(I,IQ)=-MMX(IBLOCK(I,IQ),ITYPE(I,IQ))
	MY(I,IQ)=-MMY(IBLOCK(I,IQ),ITYPE(I,IQ))
	MZ(I,IQ)=MMZ(IBLOCK(I,IQ),ITYPE(I,IQ))
C
	ENDIF
C
	IF (IOR(I,IQ).EQ.-1) THEN
	MX(I,IQ)=-MX(I,IQ)
	MZ(I,IQ)=-MZ(I,IQ)
	ENDIF
C
	ENDIF
C
	IF (ITYPE(I,IQ).EQ.2.OR.ITYPE(I,IQ).EQ.4) THEN
C
	IF (IDIR(I,IQ).EQ.1) THEN
C
	MX(I,IQ)=MMX(IBLOCK(I,IQ),ITYPE(I,IQ))
	MY(I,IQ)=MMY(IBLOCK(I,IQ),ITYPE(I,IQ))
	MZ(I,IQ)=MMZ(IBLOCK(I,IQ),ITYPE(I,IQ))
C
	ENDIF
C
	IF (IDIR(I,IQ).EQ.-1) THEN
C
	MX(I,IQ)=-MMX(IBLOCK(I,IQ),ITYPE(I,IQ))
	MY(I,IQ)=MMY(IBLOCK(I,IQ),ITYPE(I,IQ))
	MZ(I,IQ)=-MMZ(IBLOCK(I,IQ),ITYPE(I,IQ))
C
	ENDIF
C
	IF (IOR(I,IQ).EQ.-1) THEN
	MX(I,IQ)=-MX(I,IQ)
	MY(I,IQ)=-MY(I,IQ)
	ENDIF
C
	ENDIF
C
	END DO
C
	END DO
C
	RETURN
	END
C
C	-------------------------------
	SUBROUTINE PMU(C,B,I,IPOS1,IQ1)
C	-------------------------------
	REAL MX,MY,MZ,M
	COMMON/PMUDAT/X1(424,2),X2(424,2),Y1(424,2),Y2(424,2),
     $	Z1(424,2),Z2(424,2),MX(424,2),MY(424,2),MZ(424,2)
	DIMENSION C(3),V1(3),V2(3),M(3)
C
C	CALCULATE THE I'TH COMPONENT OF FIELD AT THE FIELD POINT C
C	IN AN UNDULATOR MAGNET SPECIFIED IN PMUDAT
C
	B=0.0
C
C	ADD CONTRIBUTION FROM EACH BLOCK IN TURN
C
	IQMIN=1
	IQMAX=2
	IAMIN=1
	IAMAX=NB
C
	IF (IQ1.NE.0) THEN
	IQMIN=IQ1
	IQMAX=IQ1
	ENDIF
C
	IF (IPOS1.NE.0) THEN
	IAMIN=IPOS1
	IAMAX=IPOS1
	ENDIF

	DO 40 IQ=IQMIN,IQMAX
	DO 40 IA=IAMIN,IAMAX
C
C	SET UP COORDINATES OF THE BLOCK CORNERS IN VECTORS V1,V2
C
	V1(1)=X1(IA,IQ)
	V2(1)=X2(IA,IQ)
	V1(3)=Z1(IA,IQ)
	V2(3)=Z2(IA,IQ)
	V1(2)=Y1(IA,IQ)
	V2(2)=Y2(IA,IQ)
C
C	SET UP MAGNETISATION VECTOR
C
	M(1)=MX(IA,IQ)
	M(2)=MY(IA,IQ)
	M(3)=MZ(IA,IQ)
C
C	CALCULATE FIELD FROM THIS BLOCK USING PMB ROUTINE
C
	CALL PMB_NEW(C,V1,V2,M,B,I)
40	CONTINUE
	RETURN
	END
C
C	---------------------------------
	SUBROUTINE PMB_NEW(C,V1,V2,M,B,I)
C	---------------------------------
	REAL I1,I2,M
	COMMON/PMBDAT/R1I,R1J,R1K,R2I,R2J,R2K,A1,A2,A3,A4,A5,A6,A7,A8
	DIMENSION C(3),V1(3),V2(3),M(3),R1(3),R2(3)
	DATA FOURPI/12.5663706144/
C
C	CALCULATE THE I'TH COMPONENT OF FIELD AT THE FIELD POINT C
C	DUE TO A SINGLE PERMANENT MAGNET BLOCK WITH COORDINATES V1,V2
C	AND MAGNETISATION M
C
C	SET UP RELATIVE COORDINATES
C
	DO 10 L=1,3
	R1(L)=C(L)-V1(L)
	R2(L)=C(L)-V2(L)
10	CONTINUE
C
C	EVALUATE THE INTEGRALS
C
	DO 20 J=1,3
	IF(J.EQ.I)GOTO 20
	K=6-(I+J)
	R1I=R1(I)
	R1J=R1(J)
	R2I=R2(I)
	R2J=R2(J)
	IF(R1(K).GE.0.0.AND.R2(K).GE.0.0)THEN
	  R1K=R1(K)
	  R2K=R2(K)
	ELSE
	  R1K=-R2(K)
	  R2K=-R1(K)
	ENDIF
	A1=sqrt((R2I*R2I)+(R2J*R2J)+(R2K*R2K))
	A2=sqrt((R1I*R1I)+(R1J*R1J)+(R2K*R2K))
	A3=sqrt((R1I*R1I)+(R2J*R2J)+(R1K*R1K))
	A4=sqrt((R2I*R2I)+(R1J*R1J)+(R1K*R1K))
	A5=sqrt((R1I*R1I)+(R2J*R2J)+(R2K*R2K))
	A6=sqrt((R2I*R2I)+(R1J*R1J)+(R2K*R2K))
	A7=sqrt((R2I*R2I)+(R2J*R2J)+(R1K*R1K))
	A8=sqrt((R1I*R1I)+(R1J*R1J)+(R1K*R1K))
	IF(abs(M(I)).GT.1.0e-06)B=B-(M(I)*I1()/FOURPI)
	IF(abs(M(J)).GT.1.0e-06)B=B-(M(J)*I2()/FOURPI)
20	CONTINUE
	RETURN
	END
C
C	------------------------
	FUNCTION I1()
C	------------------------
	REAL I1
	COMMON/PMBDAT/R1I,R1J,R1K,R2I,R2J,R2K,A1,A2,A3,A4,A5,A6,A7,A8

C
C	EVALUATE INTEGRAL NO. 1    R  DR  DR  / R**3
C	                            J   I   K

	I1=0.0
	IF(abs(R1J).LT.1.0D-08)GOTO 10
	B2=R1I*R2K/(R1J*A2)
	B4=R2I*R1K/(R1J*A4)
	B6=R2I*R2K/(R1J*A6)
	B8=R1I*R1K/(R1J*A8)
	I1=atan(B2)+atan(B4)-atan(B6)-atan(B8)
10	IF(abs(R2J).LT.1.0D-08)RETURN
	B1=R2I*R2K/(R2J*A1)
	B3=R1I*R1K/(R2J*A3)
	B5=R1I*R2K/(R2J*A5)
	B7=R2I*R1K/(R2J*A7)
	I1=I1+atan(B1)+atan(B3)-atan(B5)-atan(B7)
	RETURN
	END
C
C	------------------------
	FUNCTION I2()
C	------------------------
	REAL I2
	COMMON/PMBDAT/R1I,R1J,R1K,R2I,R2J,R2K,A1,A2,A3,A4,A5,A6,A7,A8

C
C	EVALUATE INTEGRAL NO. 2    R  DR  DR  / R**3
C	                            J   J   K

	B1=A1+R2K
	B2=A2+R2K
	B3=A3+R1K
	B4=A4+R1K
	B5=A5+R2K
	B6=A6+R2K
	B7=A7+R1K
	B8=A8+R1K
	I2=log((B1*B2*B3*B4)/(B5*B6*B7*B8))
	RETURN
	END
C
C	---------------------------
	SUBROUTINE METROP(DE,T,ANS)
C	---------------------------
	LOGICAL ANS
C
C	Metropolis Algorithm. ANS is a logical variable which issues a verdict
C	on whether to accept a reconfiguration which leads to a change DE in
C	the objective function E. If DE<0, ANS = .TRUE., while if DE>0, ANS is
C	only .TRUE. with probability exp(-DE/T), where T is a temperature
C	determined by the annealing schedule.
C
	IF (DE.LE.0.0) THEN
	ANS=.TRUE.
	ELSE
	IF(RAND(0).LT.EXP(-DE/T)) THEN
	ANS=.TRUE.
	ELSE
	ANS=.FALSE.
	ENDIF
	ENDIF
	RETURN
	END
C
C	----------------------
	SUBROUTINE SWAP(L1,L2)
C	----------------------
	LTEMP=L1
	L1=L2
	L2=LTEMP
	RETURN
	END
C
C	----------------
	SUBROUTINE PHASE
C	----------------
	COMMON/PAR/PERIOD,GAP,HEIGHT,WIDTH,SEP,NP,BS1,NB
	COMMON/ZPAR/ZMIN,ZMAX,DELZ,ANEX,NSTP,NZSTEP
	COMMON/OUT/FIRSTIX(-12:12),FIRSTIY(-12:12),
     $	SECNDIX(-12:12),SECNDIY(-12:12),BUNDX(0:9999),
     $	BUNDY(0:9999),PHQUAD(2),PHERR,FIRSTIXMAX,
     $	FIRSTIYMAX,SECNDIXMAX,SECNDIYMAX,STRX,STRY
	COMMON/TRAJ/X(0:9999),XP(0:9999),Y(0:9999),YP(0:9999)
	DIMENSION PH(0:9999),PH2(0:9999),V1(401),V2(401)
C
	EG=3.0
	GAMMA=EG/0.511E-3
	CONST=0.03/EG*1E-6
	TWOPI=2.0*ACOS(-1.0)
	HT=1.054588E-27
	EVTOERG=1.602E-12
	C=2.9911124E8
	TWOC=2.0*C
C
        XX1=0.0
        XX=0.0
	BA=BUNDY(0)*1E4
        X1A=0.0
        XP(0)=X1A
	X(0)=0
c
        DO I=1,NZSTEP
	BZ=BUNDY(I)*1E4
        XX1=XX1+CONST*(BA+BZ)/2.0*DELZ
        X1Z=XX1
        XP(I)=XX1
	X(I)=X(I-1)+(XP(I)+XP(I-1))/2*DELZ
        BA=BZ
        X1A=X1Z
        END DO
C
        YY1=0.0
        YY=0.0
	BA=BUNDX(0)*1E4
        Y1A=0.0
	YP(0)=Y1A
	Y(0)=0.0
C
        DO I=1,NZSTEP
	BZ=BUNDX(I)*1E4
        YY1=YY1-CONST*(BA+BZ)/2.0*DELZ
        Y1Z=YY1
        YP(I)=YY1
	Y(I)=Y(I-1)+(YP(I)+YP(I-1))/2*DELZ
        BA=BZ
        Y1A=Y1Z
        END DO
C
	THETAX=0.0
	THETAY=0.0
C
	WX=0.0
        WXA=0.0
	WY=0.0
        WYA=0.0
	PH(0)=0.0
	DO I=1,NZSTEP
	WXZ=(XP(I)-THETAX)**2
	WYZ=(YP(I)-THETAY)**2
	WX=WX+(WXA+WXZ)/2*DELZ*1E-3
	WY=WY+(WYA+WYZ)/2*DELZ*1E-3
	PH(I)=(WX+WY)/TWOC
	PH2(I)=(I*DELZ*1E-3/GAMMA/GAMMA/TWOC+PH(I))
	WXA=WXZ
	WYA=WYZ
        END DO
C
C       Compute fundamental energy
C
	NPOLMAX=4*NP+1
	NSKIP=8
	J1=NSKIP+1
	J2=NPOLMAX-NSKIP
	NDATA=NPOLMAX-2*NSKIP
C
	DO J=J1,J2
	JJ=J-NSKIP
	IZP=NZSTEP/2-NP*NSTP/2+(J-1)*NSTP/4
	V1(JJ)=IZP
	V2(JJ)=PH2(IZP)
	END DO
	CALL FIT(V1,V2,NDATA,SIG,0,A,B,SIGA,SIGB,CHI2)
	OMEGA0=TWOPI/B/NSTP
C	WRITE (6,*)'OMEGA0 =',OMEGA0
C	WRITE (6,*)'E1 (eV)=',OMEGA0*HT/EVTOERG
C
C       Compute rms phase error
C
	DO J=J1,J2
	JJ=J-NSKIP
	IZP=NZSTEP/2-NP*NSTP/2+(J-1)*NSTP/4
	V1(JJ)=IZP
	V2(JJ)=PH(IZP)
	END DO
	CALL FIT(V1,V2,NDATA,SIG,0,A,B,SIGA,SIGB,CHI2)
C
	DO J=J1,J2
	IZP=NZSTEP/2-NP*NSTP/2+(J-1)*NSTP/4
	PHFIT=A+B*IZP
	PH(IZP)=PH(IZP)-PHFIT
C	Z=ZMIN+(ZMAX-ZMIN)*IZP/NZSTEP
C	WRITE(7,*)Z,PH(IZP)*OMEGA0*360/TWOPI
	END DO
C
	PHERR=0.0
	DO J=J1,J2
	IZP=NZSTEP/2-NP*NSTP/2+(J-1)*NSTP/4
	PHERR=PHERR+PH(IZP)*OMEGA0*PH(IZP)*OMEGA0
	END DO
	PHERR=SQRT(PHERR/NDATA)*360/TWOPI
C
	RETURN
	END
C
C	-------------------
	SUBROUTINE STRAIGHT
C	-------------------
	COMMON/PAR/PERIOD,GAP,HEIGHT,WIDTH,SEP,NP,BS1,NB
	COMMON/ZPAR/ZMIN,ZMAX,DELZ,ANEX,NSTP,NZSTEP
	COMMON/OUT/FIRSTIX(-12:12),FIRSTIY(-12:12),
     $	SECNDIX(-12:12),SECNDIY(-12:12),BUNDX(0:9999),
     $	BUNDY(0:9999),PHQUAD(2),PHERR,FIRSTIXMAX,
     $	FIRSTIYMAX,SECNDIXMAX,SECNDIYMAX,STRX,STRY
	COMMON/TRAJ/X(0:9999),XP(0:9999),Y(0:9999),YP(0:9999)
C
	NPOLMAX=2*NP+1
	NSKIP=8
	J1=NSKIP/2+1
	J2=NPOLMAX-NSKIP/2
	NCOUNT=J2-J1+1
C
	STRX=0.0
	STRY=0.0
	XMEAN=0.0
	DO J=J1,J2
	IZP=NZSTEP/2-(NP*NSTP)/2+(J-1)*NSTP/2+NSTP/4
	XMEAN=XMEAN+X(IZP)
	END DO
	XMEAN=XMEAN/NCOUNT
	DO J=J1,J2
	IZP=NZSTEP/2-(NP*NSTP)/2+(J-1)*NSTP/2+NSTP/4
	DXMEAN=ABS(X(IZP)-XMEAN)
	IF (DXMEAN.GT.STRX) STRX=DXMEAN
	IF (ABS(Y(IZP)).GT.STRY) STRY=ABS(Y(IZP))
	END DO
C
	RETURN
	END
C
C     ----------------------------------------------------
      SUBROUTINE FIT(X,Y,NDATA,SIG,MWT,A,B,SIGA,SIGB,CHI2)
C     ----------------------------------------------------
      DIMENSION X(NDATA),Y(NDATA),SIG(NDATA)
      SX=0.
      SY=0.
      ST2=0.
      B=0.
      IF(MWT.NE.0) THEN
        SS=0.
        DO 11 I=1,NDATA
          WT=1./(SIG(I)**2)
          SS=SS+WT
          SX=SX+X(I)*WT
          SY=SY+Y(I)*WT
11      CONTINUE
      ELSE
        DO 12 I=1,NDATA
          SX=SX+X(I)
          SY=SY+Y(I)
12      CONTINUE
        SS=FLOAT(NDATA)
      ENDIF
      SXOSS=SX/SS
      IF(MWT.NE.0) THEN
        DO 13 I=1,NDATA
          T=(X(I)-SXOSS)/SIG(I)
          ST2=ST2+T*T
          B=B+T*Y(I)/SIG(I)
13      CONTINUE
      ELSE
        DO 14 I=1,NDATA
          T=X(I)-SXOSS
          ST2=ST2+T*T
          B=B+T*Y(I)
14      CONTINUE
      ENDIF
      B=B/ST2
      A=(SY-SX*B)/SS
      SIGA=SQRT((1.+SX*SX/(SS*ST2))/SS)
      SIGB=SQRT(1./ST2)
      CHI2=0.
      IF(MWT.EQ.0) THEN
        DO 15 I=1,NDATA
          CHI2=CHI2+(Y(I)-A-B*X(I))**2
15      CONTINUE
        Q=1.
        SIGDAT=SQRT(CHI2/(NDATA-2))
        SIGA=SIGA*SIGDAT
        SIGB=SIGB*SIGDAT
      ELSE
        DO 16 I=1,NDATA
          CHI2=CHI2+((Y(I)-A-B*X(I))/SIG(I))**2
16      CONTINUE
      ENDIF
      RETURN
      END
C
C	-----------------------
	SUBROUTINE MOVE1(ISEED)
C	-----------------------
	COMMON/PAR/PERIOD,GAP,HEIGHT,WIDTH,SEP,NP,BS1,NB
	COMMON/CONF/ITYPE(424,2),IDIR(424,2),IOR(424,2),IBLOCK(424,2),
     $	IUSE(999,4),JBLOCK(999,4),NAEXT,NBEXT,NCEXT,NDEXT
C
C	SWAP TWO A-BLOCKS
C
	IF (ISEED.EQ.0) GOTO 33
C
	IMIN=4
	IMAX=(NB-4)/2
C
	N1=INT(RAND(0)*(IMAX-IMIN+1.0)/2)*2+IMIN
	IQ1=INT(RAND(0)*2.0+1.0)
1	N2=INT(RAND(0)*(IMAX-IMIN+1.0)/2)*2+IMIN
	IQ2=INT(RAND(0)*2.0+1.0)
	IF (N1.EQ.N2.AND.IQ2.EQ.IQ1) GOTO 1
C
	IF (N2.EQ.N1.AND.IQ1.EQ.IQ2) THEN
	WRITE(6,*)'*** ERROR IN MOVE1 ***'
	WRITE(6,*)'SWAPPING BLOCK ',IBLOCK(N1,IQ1)
	WRITE(6,*)'POSITION ',N1,' ARRAY ',IQ1,' TYPE ',ITYPE(N1,IQ1)
	WRITE(6,*)'WITH BLOCK ',IBLOCK(N2,IQ2)
	WRITE(6,*)'POSITION ',N2,' ARRAY ',IQ2,' TYPE ',ITYPE(N2,IQ2)
	STOP
	ENDIF
C
33	CALL SWAP(IBLOCK(N1,IQ1),IBLOCK(N2,IQ2))
C
	RETURN
	END
C
C	-----------------------
	SUBROUTINE MOVE2(ISEED)
C	-----------------------
	COMMON/PAR/PERIOD,GAP,HEIGHT,WIDTH,SEP,NP,BS1,NB
	COMMON/CONF/ITYPE(424,2),IDIR(424,2),IOR(424,2),IBLOCK(424,2),
     $	IUSE(999,4),JBLOCK(999,4),NAEXT,NBEXT,NCEXT,NDEXT
C
C	SWAP TWO B-BLOCKS
C
	IF (ISEED.EQ.0) GOTO 33
C
	IMIN=3
	IMAX=NB-2
C
	N1=INT(RAND(0)*(IMAX-IMIN+1.0)/2)*2+IMIN
	IQ1=INT(RAND(0)*2.0+1.0)
1	N2=INT(RAND(0)*(IMAX-IMIN+1.0)/2)*2+IMIN
	IQ2=INT(RAND(0)*2.0+1.0)
	IF (N2.EQ.N1.AND.IQ1.EQ.IQ2) GOTO 1
C
	IF (N2.EQ.N1.AND.IQ1.EQ.IQ2) THEN
	WRITE(6,*)'*** ERROR IN MOVE2 ***'
	WRITE(6,*)'SWAPPING BLOCK ',IBLOCK(N1,IQ1)
	WRITE(6,*)'POSITION ',N1,' ARRAY ',IQ1,' TYPE ',ITYPE(N1,IQ1)
	WRITE(6,*)'WITH BLOCK ',IBLOCK(N2,IQ2)
	WRITE(6,*)'POSITION ',N2,' ARRAY ',IQ2,' TYPE ',ITYPE(N2,IQ2)
	STOP
	ENDIF
C
33	CALL SWAP(IBLOCK(N1,IQ1),IBLOCK(N2,IQ2))
C
	RETURN
	END
C
C	-----------------------
	SUBROUTINE MOVE3(ISEED)
C	-----------------------
	COMMON/PAR/PERIOD,GAP,HEIGHT,WIDTH,SEP,NP,BS1,NB
	COMMON/CONF/ITYPE(424,2),IDIR(424,2),IOR(424,2),IBLOCK(424,2),
     $	IUSE(999,4),JBLOCK(999,4),NAEXT,NBEXT,NCEXT,NDEXT
C
C	SWAP TWO C-BLOCKS
C
	IF (ISEED.EQ.0) GOTO 33
C
	NN1=INT(RAND(0)*2.0)
	IQ1=INT(RAND(0)*2.0+1.0)
1	NN2=INT(RAND(0)*2.0)
	IQ2=INT(RAND(0)*2.0+1.0)
	IF (NN2.EQ.NN1.AND.IQ1.EQ.IQ2) GOTO 1
c
	IF (NN1.EQ.0) N1=2
	IF (NN1.EQ.1) N1=NB-1
	IF (NN2.EQ.0) N2=2
	IF (NN2.EQ.1) N2=NB-1
C
	IF (N2.EQ.N1.AND.IQ1.EQ.IQ2) THEN
	WRITE(6,*)'*** ERROR IN MOVE3 ***'
	WRITE(6,*)'SWAPPING BLOCK ',IBLOCK(N1,IQ1)
	WRITE(6,*)'POSITION ',N1,' ARRAY ',IQ1,' TYPE ',ITYPE(N1,IQ1)
	WRITE(6,*)'WITH BLOCK ',IBLOCK(N2,IQ2)
	WRITE(6,*)'POSITION ',N2,' ARRAY ',IQ2,' TYPE ',ITYPE(N2,IQ2)
	STOP
	ENDIF
C
33	CALL SWAP(IBLOCK(N1,IQ1),IBLOCK(N2,IQ2))
C
	RETURN
	END
C
C	-----------------------
	SUBROUTINE MOVE4(ISEED)
C	-----------------------
	COMMON/PAR/PERIOD,GAP,HEIGHT,WIDTH,SEP,NP,BS1,NB
	COMMON/CONF/ITYPE(424,2),IDIR(424,2),IOR(424,2),IBLOCK(424,2),
     $	IUSE(999,4),JBLOCK(999,4),NAEXT,NBEXT,NCEXT,NDEXT
C
C	SWAP TWO D-BLOCKS
C
	IF (ISEED.EQ.0) GOTO 33
C
	NN1=INT(RAND(0)*2.0)
	IQ1=INT(RAND(0)*2.0+1.0)
1	NN2=INT(RAND(0)*2.0)
	IQ2=INT(RAND(0)*2.0+1.0)
	IF (NN2.EQ.NN1.AND.IQ1.EQ.IQ2) GOTO 1
C
	IF (NN1.EQ.0) N1=1
	IF (NN1.EQ.1) N1=NB
	IF (NN2.EQ.0) N2=1
	IF (NN2.EQ.1) N2=NB
C
	IF (N2.EQ.N1.AND.IQ1.EQ.IQ2) THEN
	WRITE(6,*)'*** ERROR IN MOVE4 ***'
	WRITE(6,*)'SWAPPING BLOCK ',IBLOCK(N1,IQ1)
	WRITE(6,*)'POSITION ',N1,' ARRAY ',IQ1,' TYPE ',ITYPE(N1,IQ1)
	WRITE(6,*)'WITH BLOCK ',IBLOCK(N2,IQ2)
	WRITE(6,*)'POSITION ',N2,' ARRAY ',IQ2,' TYPE ',ITYPE(N2,IQ2)
	STOP
	ENDIF
C
33	CALL SWAP(IBLOCK(N1,IQ1),IBLOCK(N2,IQ2))
C
	RETURN
	END
C
C	-----------------------
	SUBROUTINE MOVE5(ISEED)
C	-----------------------
	COMMON/PAR/PERIOD,GAP,HEIGHT,WIDTH,SEP,NP,BS1,NB
	COMMON/CONF/ITYPE(424,2),IDIR(424,2),IOR(424,2),IBLOCK(424,2),
     $	IUSE(999,4),JBLOCK(999,4),NAEXT,NBEXT,NCEXT,NDEXT
C
C	REVERSE ORIENTATION OF A BLOCK (ANY TYPE)
C
	IF (ISEED.EQ.0) GOTO 33
C
	N=INT(RAND(0)*NB+1.0)
	IQ=INT(RAND(0)*2.0+1.0)
C
33	IOR(N,IQ)=-IOR(N,IQ)
C
	RETURN
	END
C
C	-----------------------
	SUBROUTINE MOVE6(ISEED)
C	-----------------------
	COMMON/PAR/PERIOD,GAP,HEIGHT,WIDTH,SEP,NP,BS1,NB
	COMMON/CONF/ITYPE(424,2),IDIR(424,2),IOR(424,2),IBLOCK(424,2),
     $	IUSE(999,4),JBLOCK(999,4),NAEXT,NBEXT,NCEXT,NDEXT
C
C	REPLACE A RANDOM BLOCK WITH A SPARE ONE
C
	IF (ISEED.EQ.0) GOTO 33
C
	N=INT(RAND(0)*NB+1.0)
	IQ=INT(RAND(0)*2.0+1.0)
	IT=ITYPE(N,IQ)
C
	IF (IT.EQ.1) NNMAX=NAEXT
	IF (IT.EQ.2) NNMAX=NBEXT
	IF (IT.EQ.3) NNMAX=NCEXT
	IF (IT.EQ.4) NNMAX=NDEXT
C
	NN=INT(RAND(0)*NNMAX+1.0)
C
33	CALL SWAP(IBLOCK(N,IQ),JBLOCK(NN,IT))
C
	RETURN
	END

ext: for

max_stars_repo_name: DiamondLightSource/Opt-ID

max_issues_repo_name: DiamondLightSource/Opt-ID

max_forks_repo_name: JossWhittle/Opt-ID


text:
      subroutine zone_nnum
!***********************************************************************
! $Id: zone_nnum.f,v 1.2 2006/08/16 21:37:29 zvd Exp $
!***********************************************************************
!  Copyright, 2002, 2004,  The  Regents of the University of California.
!  This  program  was  prepared by  the  Regents  of the  University  of
!  California  at Los Alamos National Laboratory (the  University) under
!  contract No. W-7405-ENG-36 with the U. S. Department of Energy (DOE).
!  All rights in the program are reserved by the DOE and the University.
!  Permission is  granted to the  public to copy  and  use this software
!  without  charge,  provided that this  Notice  and  any  statement  of
!  authorship are reproduced on all copies. Neither the U. S. Government
!  nor  the  University  makes any  warranty,  express  or  implied,  or
!  assumes any liability or responsibility for the use of this software.
!***********************************************************************

      use comunits, only : rock_unit_number
      use comgrid, only : n_grid_points, izonef
      implicit none

      character*80 wdd1

      logical null1
      integer cnum, i, izone, izonel, nin, nodez, nsl
      integer nxy,icnl_old, zone_unit_number, open_file
      character* 4 macro
      character*80 ltest
      character*100 zone_file_name
      real*8 xg, xz(8), yg, yz(8), zg, zz(8)
      real*8 tol_zone, zxy_min, zxy_max
      integer imodel, j, n_n_n
      integer zone_dpadd

      integer, allocatable :: ncord(:)

      allocate(ncord(n_grid_points))

      read (rock_unit_number, '(a4)', end=90)  macro
      if (macro .eq. 'file') then
! zone definitions are in a separate file
         read (rock_unit_number, '(a100)') zone_file_name
         zone_unit_number = open_file(zone_file_name, 1, 'OLD')
         do
            read (zone_unit_number, '(a4)', end=90) macro
            if (macro .eq. 'zone') exit
         end do
      else
         zone_unit_number = rock_unit_number
         backspace zone_unit_number
      end if

      do
         read (zone_unit_number, '(a80)', end=90) wdd1
         if (null1(wdd1) .or. wdd1(1:4) .eq. 'rock') exit
         read(wdd1, *) izone
         read(zone_unit_number, *) macro
         if(macro .eq. 'list') then
! zone input using list format
! read in coordinates for nodes in zone
            nin = 0
            do
               read(zone_unit_number, '(a80)', end = 90) ltest
               if(.not.null1(ltest)) exit
               read(ltest, *, end = 80, err = 80) xg, yg, zg
               nin = nin + 1
               call near3(xg, yg, zg, nodez)
               ncord(nin) = nodez
               izonef(nodez) = izone
            end do
 80         continue
         else if (macro .eq. 'nnum') then
! zone input using nnum format
            read (zone_unit_number, *) nin, (ncord (i), i = 1, nin)
            do i = 1, nin
               izonef(ncord(i)) = izone
            end do
         endif
      end do
 90   continue
      deallocate(ncord)
      if (zone_unit_number .ne. rock_unit_number)
     &     close(zone_unit_number)
      return
      end

max_stars_repo_name: lanl/PLUMECALC

max_issues_repo_name: lanl/PLUMECALC

max_forks_repo_name: lanl/PLUMECALC


text:
!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------

!C
!C***
!C*** module hecmw_precond_ML_nn
!C***
!C
module hecmw_precond_ML_nn
  use  hecmw_util

  private

  public:: hecmw_precond_ML_nn_setup
  public:: hecmw_precond_ML_nn_apply
  public:: hecmw_precond_ML_nn_clear

  integer(kind=kint), save :: id

  logical, save :: INITIALIZED = .false.

contains

  subroutine hecmw_precond_ML_nn_setup(hecMAT, hecMESH, sym)
    use hecmw_matrix_misc
    use hecmw_mat_id
    implicit none
    type(hecmwST_matrix), intent(inout) :: hecMAT
    type(hecmwST_local_mesh), intent(in) :: hecMESH
    integer(kind=kint), intent(in) :: sym
    integer(kind=kint) :: ierr
    integer(kind=kint), save :: n_recycle = 0
    if (INITIALIZED) then
      if (hecMAT%Iarray(98) == 1) then      ! need symbolic and numerical setup
        call hecmw_precond_ML_nn_clear()
      else if (hecMAT%Iarray(97) == 1) then ! need numerical setup only
        call hecmw_precond_ML_nn_clear()
      else                                  ! no need to setup or skip setup
        call hecmw_mat_id_clear(id)
        call hecmw_mat_id_set(hecMAT, hecMESH, id)
        return
      endif
    endif
    call hecmw_mat_id_set(hecMAT, hecMESH, id)
    call hecmw_ML_wrapper_setup(id, sym, hecMAT%NDOF, ierr)
    INITIALIZED = .true.
    hecMAT%Iarray(98) = 0 ! symbolic setup done
    hecMAT%Iarray(97) = 0 ! numerical setup done
    n_recycle = 0
  end subroutine hecmw_precond_ML_nn_setup

  subroutine hecmw_precond_ML_nn_apply(WW)
    implicit none
    real(kind=kreal), intent(inout) :: WW(:)
    integer(kind=kint) :: ierr
    call hecmw_ML_wrapper_apply(id, WW, ierr)
  end subroutine hecmw_precond_ML_nn_apply

  subroutine hecmw_precond_ML_nn_clear()
    use hecmw_mat_id
    implicit none
    integer(kind=kint) :: ierr
    call hecmw_ML_wrapper_clear(id, ierr)
    call hecmw_mat_id_clear(id)
    INITIALIZED = .false.
  end subroutine hecmw_precond_ML_nn_clear

end module     hecmw_precond_ML_nn

max_stars_repo_name: masae-hayashi/FrontISTR

max_issues_repo_name: masae-hayashi/FrontISTR

max_forks_repo_name: masae-hayashi/FrontISTR