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
