	SUBROUTINE CENTER(IFND,ityp,width,tthzero)
C
C ROUTINE TO CENTER A REFLECTION IN THE DETECTOR.  REFLECTION MUST BE
C  CLOSE BEFORE CALLING THIS ROUTINE.  IF MAXIMUM COUNT RATE IS BELOW
C  50 CPS, IFND WILL BE NON-ZERO UPON RETURN: also set if ierr flag set by drive
C ITYPE IS NON-ZERO WHEN 8-REFLECTION MODE IS IN EFFECT
c WIDTH is the GAMMA for the reflection
C
C  NOTE: THIS ROUTINE WILL LEAVE THE SHUTTER OPEN
C
	include 'single.inc'
C
c 20-jul-1993 added wait statement 
c
	call wait
	CALL MTSTEP(2,0.10,CMAX,CMIN,.5,ns,ierr)	!INITIAL OMEGA SCAN
c send failure indication if peak not found
	if(cmax.lt.minct)go to 1000
	if(ierr.eq.1)goto1000
c also fail if high background but no peak (cmax/cmin < 4)
	if(cmin/cmax.gt.1./maxmin)go to 1000
c
c do parabola scan to check if peak not streak and center 2theta
c
	CALL MTSTEP(1,.1,CMAX,CMIN,.5,ns,ierr)	!2THETA SCAN
	if(ierr.eq.1)goto1000
	if(cmax .lt. minct)goto1000
	if(cmin/cmax.gt. 1./maxmin)goto1000
c position omega at peak
	call mtstep(2,0.05,cmax,cmin,.5,ns,ierr)
	if(ierr.eq.1)goto1000
	if(cmax.lt.minct)go to 1000
	call mtstep(3,.3,cmax,cmin,.5,ns,ierr)	!Initial Chi scan
	if(ierr.eq.1)goto1000
	if(cmax.lt.minct)go to 1000
	call minicenter(ifnd,width)		!Center chi and omega
	if(ifnd.ne.0)go to 1000
	if(ityp.ne.0)return
c Prepare to measure on -2theta side
	call ofbeam
	omega=rdmtr(2)
	tth=rdmtr(1)
	thpos=omega+0.5*tth
	call drive(1,-tth,0,1,ierr)
	if(ierr .eq. 1)goto1000
	call wait
	call onbeam
	CALL MTSTEP(2,0.10,CMAX,CMIN,.5,ns,ierr)	!INITIAL OMEGA SCAN
	if(ierr.eq.1)goto1000
	if(cmax.lt.minct)go to 1000
	T=555.5/CMAX
	IF(T.GT.tmax)T=tmax
	if(t.lt.0.25)t=0.25
	CALL MTSTEP(1,.1,CMAX,CMIN,t,ns,ierr)	!2THETA SCAN
	if(cmax.lt.minct)go to 1000
	T=555.5/CMAX
	IF(T.GT.tmax)T=tmax
	if(t.lt.0.25)t=0.25
	call omstep(ifnd,pos,t,xwid)
	if(ifnd.ne.0)return
	call ofbeam
c
c calc final position of peak
c
	thminus=pos+0.5*rdmtr(1)
	tth=thpos-thminus
	omega=(pos+omega)*0.5
        tthzero=thpos+thminus-2.*omega
c drive two-theta to maximum, allowing for zero error
	call drive(1,tth+tthzero,0,1,ierr)
	if(ierr .eq. 1)goto1000
	call drive(2,omega,0,1,ierr)
	if(ierr .eq. 1)goto1000
	call wait
	call onbeam
	return
c Weak Peak or Error - Set Flag
1000	ifnd=1
	call ofbeam
	return
	end
	subroutine minicenter(ifnd,width)
c
c routine to do omega parabolic centering, chi edge finding, and
c   omega profiling for a reflection
c
c 'ifnd' is non-zero if reflection not found 
c       (also set if out-of-bounds flag set by drive)
c 'width' is reflection width (in omega)
c
c do omega parabolic scan
c
	include 'single.inc'
	CALL MTSTEP(2,0.05,CMAX,CMIN,.5,ns,ierr)
	if(ierr.eq.1)goto1000
	if(cmax.lt.minct)go to 1000
	T=555.5/CMAX
	IF(T.GT.tmax)T=tmax
	if(t.lt.0.25)t=0.25
	PI=RDMTR(3)
	CMAX=0.75*CMAX
	DEL=0.32
C GET POSITIVE EDGE OF REFLECTION FOR CHI
	start=pi
	CALL CONV(3,PI,DEL,CMAX,T,ierr)
	if(ierr .eq. 1)goto1000
        call drive(3,start,0.0,0,ierr)
	TEMP=PI
	DEL=-0.32
C NOW GET MINUS CHI EDGE
	CALL CONV(3,PI,DEL,CMAX,T,ierr)
	if(ierr .eq. 1)goto1000
	TEMP=(TEMP+PI)*0.5
	CALL DRIVE(3,TEMP,0.0,1,ierr)
	if(ierr .eq. 1)goto1000
C PERFORM THETA STEP SCAN: BRANCH ON LOGICAL BADCRYSTAl
	if(badcrystal)then
	  call mtstep(2,0.02,CMAX,CMIN,.5,ns,ierr)
	  if(ierr .ne. 0)return
	  pos=rdmtr(2)
	  width=0.0
	else
	  call omstep(ifnd,pos,t,width)
	  if(ifnd.ne.0)return
C DRIVE TO FINAL POSITION ON OMEGA
	  CALL DRIVE(2,POS,0,1,ierr)
	  if(ierr .eq. 1)goto1000
	  call wait
	  ifnd=0
	endif
	return
c here if some kind of error
1000	ifnd=1
	return
	end
	subroutine omstep(ifnd,pos,t,width)
c
c Perform 25 step Omega Scan
c
	implicit none
	real*4 PANGL(25),PCOUNT(25),PARAM(4),pos,t,width
     1,rdmtr,cmax,stpsiz,omgpos,rdsclr,rdtime,tth
	integer*4 ifnd,i,ierr
c
	OMGPOS=RDMTR(2)
	TTH=RDMTR(1)
	IF(T.GT.1.0)T=1.0
	CMAX=0.0
c Calculate step size for omega step scan
	stpsiz=0.075*abs(tan(tth/114.592))
	stpsiz=anint(100.0*stpsiz)*0.01
	if(stpsiz.lt.0.02)stpsiz=0.02
	DO 100 I=1,25
	  PANGL(I)=(I-13)*stpsiz+OMGPOS
	  CALL DRIVE(2,PANGL(I),0.0,0,ierr)
	  CALL WAIT
	  CALL COUNT(T)
	  CALL sclwt
	  PCOUNT(I)=RDSCLR()/RDTIME()
	  IF(PCOUNT(I).LE.CMAX)GO TO 100
	    CMAX=PCOUNT(I)
	    PARAM(4)=PANGL(I)
100	CONTINUE
	PARAM(1)=PCOUNT(1)
	PARAM(2)=3.*stpsiz
	PARAM(3)=CMAX*PARAM(2)/stpsiz
	CALL PKLSQ(PANGL,PCOUNT,25,PARAM,IFND)
	IF(IFND.NE.0)RETURN
	pos=param(4)
	width=param(2)
	return
	end
	SUBROUTINE CONV(n,PI,D,C,T,ierr)
C ROUTINE TO DETERMINE EDGE OF REFLECTION FOR MOTOR N
	implicit none
	include 'single.inc'
	real*4 pi,c,d,t,rdmtr,rdsclr,rdtime,dmin,start,sgn
     1,	omg,cp
	integer*4 n,ierr,ind
C
	IND=0
	ierr=0
	DMIN=0.009
	SGN=1.0
	start=pi
C ADJUST position
100	PI=PI+SGN*D
	if(n.eq.3)then
c move chi
	  CALL DRIVE(3,PI,0.0,0,ierr)
	else
c move 2theta with crystal fixed by compensating move on omega
	  omg=rdmtr(2)-0.5*sgn*d
	  call drive(2,omg,0.0,0,ierr)
	  call drive(1,pi,0.0,0,ierr)
	  if(abs(pi-start).gt.2.5)then
		ierr=1
		return
	  endif
	endif
120	CALL WAIT
	CALL COUNT(T)
	CALL sclwt
	CP=RDSCLR()/RDTIME()
	IF(CP.LT.C)GO TO 140
	SGN=1.0
	IF(IND)150,100,150
140	SGN=-1.0
	IND=1
150	D=D*0.5
	IF(ABS(D).GT.DMIN)GO TO 100
	RETURN
	END
	SUBROUTINE PROFST
C
C DOES PEAK PROFILE AT CURRENT POSITION
C
	real*4 tangl(4),TH(3)
	character*6 angname(4)
	data angname/'TH/2TH',' OMEGA','   CHI','   PHI'/
c
C FIRST WAIT FOR DRIVE COMPLETION AND CHECK NOT AT ZERO
c
   	call wait
   	call rdallmtr(tangl,TH)
	IF(ABS(tangl(1)).LT. 1.)THEN
	  TYPE *,' PROFILING PROHIBITED CLOSE TO DIRECT BEAM'
	  RETURN
	ENDIF
C
C REQUEST PROFILING PARAMETERS FROM USER
C
	TYPE *,' PEAK PROFILING'
7	TYPE 8
8	FORMAT('$INPUT 1 FOR THETA/2THETA, 2 FOR OMEGA,',
     1  ' 3 FOR CHI, 4 FOR PHI: ')
	READ(5,12,END=9999,ERR=9)IPROF
12	FORMAT(I1)
	IF(IPROF .EQ. 0)RETURN
	IF(IPROF.GT. 0 .AND. IPROF .LT. 5)GOTO19
9	TYPE *,'  INVALID VALUE: '
	GOTO7
C
19	TYPE 15
15      FORMAT('$ENTER NUMBER OF STEPS:')
	read(5,20,end=9999,err=19)NSTEP
20	FORMAT(I10)
	IF(NSTEP .EQ. 0)RETURN
	IF(NSTEP .GT. 5 .AND. NSTEP .LE. 99)GOTO29
        TYPE 25,NSTEP
25	FORMAT(' *****NSTEP VALUE OF ',I5,' NOT ALLOWED')
	GOTO19
C
29	TYPE 30
30	FORMAT('$ENTER STEP SIZE (DEGREES):')
	read(5,40,end=9999,err=29)SSIZE
40	FORMAT(F10.0)
 	IF(ABS(SSIZE) .GE. 0.005 .AND. ABS(SSIZE) .LT. 1.0)GOTO59
	TYPE 50,SSIZE
50	FORMAT(' *****STEP SIZE VALUE OF ',F10.3,' NOT ALLOWED')
	GOTO29
C
59	TYPE 60
60	FORMAT('$ENTER COUNT TIME PER STEP:')
	read(5,40,end=9999,err=59)TIME
	IF(TIME .GT. 0.05 .AND. TIME .LT. 99.)GOTO100
	TYPE 70,TIME
70	FORMAT(' *****TIME PER STEP VALUE OF ',F10.3,' SECS NOT ALLOWED')
	GOTO59
C
C PRINT HEADER FOR PROFILING
C
100	type 110,angname(iprof),tangl
110	format(/7x,a6,' PROFILE ABOUT ',4F8.2)
c	
        CALL HCALC(TANGL,TH)
	I=ABS(TH(1))+ABS(TH(2))+ABS(TH(3))
	IF(I  .GT. 0.1)THEN
   	   TYPE 120,TH
120	   FORMAT(7X,'INDICES: ',3F6.2)
	ENDIF
	CALL PKPROF(NSTEP,SSIZE,TIME,IPROF,ierr)
	if(ierr .eq. 0)return
	type 130
130	format(' Reflection OUT-OF-BOUNDS or Drive Error')
9999	RETURN
	END
C
C	 
	SUBROUTINE PKPROF(NSTEP,SSIZE,TIME,IPROF,ierr)
	parameter (maxstp=200)
C
C PERFORMS OMEGA OR THETA/2THETA STEP SCAN AND PRINTS PROFILE ABOUT CURRENT
C POSITION
C ********CHI and PHI SCANS ADDED RJA 6-NOV-93: IPROF=3,4
C
C PROFILE PARAMETERS:
C	NSTEP - NO. OF STEPS IN SCAN
C	TIME  - TIME PER STEP
C	SSIZE - STEP SIZE
C  	IPROF= 1  FOR THETA/TWO-THETA,2 FOR OMEGA. 3 FOR CHI 4 for phi
C
	include 'single.inc'
	REAL PANGL(maxstp),PCOUNT(maxstp)
	character*6 angname(4)
	data angname/'2THETA',' OMEGA','   CHI','   PHI'/
c
c Check that number of steps is possible
c
	if(nstep.gt.maxstp)then
		type *,'Too many steps in PKPROF'
		ierr=1
		return
	endif
C
C DO SCAN: NOTE THAT VARIABLES NAMED IN TERMS OF OMEGA
C
	ierr=0
	OMCENT=RDMTR(IPROF)
	CALL ONBEAM
	CMAX=0.
	IMSTP=(NSTEP+1)/2
	DO 200 I=1,NSTEP
		PANGL(I)=(I-IMSTP)*SSIZE+OMCENT
		CALL DRIVE(IPROF,PANGL(I),0,0,ierr)
		CALL WAIT
		CALL COUNT(TIME)
		CALL sclwt
		PCOUNT(I)=RDSCLR()/RDTIME()
                IF(PCOUNT(I) .GT. CMAX)THEN
			CMAX=PCOUNT(I)
			PMAX=PANGL(I)
		ENDIF
200	CONTINUE
	CALL DRIVE(IPROF,OMCENT,0,0,ierr)
	CALL OFBEAM
C
C PRINT PROFILE HEADER
C
	NTSTEP=max((NSTEP-1)/4,1)
	CCOL=CMAX/71
	TYPE 251,NINT(CMAX),ANGNAME(IPROF),PMAX
251	FORMAT(//7x,'PEAK PROFILE',15('.'),'MAX COUNT RATE = ',I6,
     1                ' CPS AT ',A6,' = ',F6.2/7X,73(1H*))
C
C PRINT PROFILE
C
	DO 300 I=1,NSTEP
		ICOL=(PCOUNT(I)/CCOL)-1
		IF(ICOL .LT. 1)ICOL=0
		ITMP=MOD(I,NTSTEP)
		IF(ITMP .EQ. 1)THEN
			TYPE 260,PANGL(I)
260			FORMAT(1X,F6.2,'|',<ICOL>X,'*',<70-ICOL>X,'|')
		ELSE
			TYPE 270
270			FORMAT(7X,'|',<ICOL>X,'*',<70-ICOL>X,'|')
		ENDIF
300	CONTINUE
	TYPE 280
280	FORMAT(7X,73(1H*)//)
C
C
	RETURN
        END
	subroutine mtstep(no,ddelo,rtmax,rtmin,time,nstep,ierr)
c
c Routine to perform a step scan of 9 steps.
c
c 'NO' is motor number (1 = theta-2theta scan)
c 'DDELO' is delta angle between steps
c 'RTMAX' is the maximum count rate
c 'RTMIN' is the minimum count rate
c 'TIME' is count time in seconds.
c 'NSTEP' is Number of Step at which peak occurs
c
	include 'single.inc'
c
c Make sure driving done
c
	ierr=0
	call wait
c drive to initial position
	delo=ddelo				!need local copy of delo
	posmax=rdmtr(no)
	ith=0
c	if(no.eq.1)then
c	  ith=1
c	  omeg0=rdmtr(2)
c	  tth0=rdmtr(1)
c	endif
	nloop=0
90	pos=posmax-4*delo
	curpos=posmax
	rtmin=1.0e6
	call drive(no,pos,0,1,ierr)
	  if(ierr .ne. 0)return
c	if(ith.ne.0)call drive(2,omeg0-(pos-tth0)*0.5,0,1,ierr)
	  if(ierr .ne. 0)return
	call wait
	call onbeam
	rtmax=0.0
	sumy=0.0
	sumxy=0.0
	sumx2y=0.0
	do 110 i=1,9
	  call count(time)
	  call sclwt
	  rate=rdsclr()/rdtime()
	  peak=rate
	  sumy=sumy+peak
	  sumxy=sumxy+(i-5)*peak
	  sumx2y=sumx2y+(i-5)**2*peak
	  if(rate.lt.rtmin)rtmin=rate
	  if(rate.le.rtmax)go to 100
	    rtmax=rate
	    posmax=rdmtr(no)
	    nstep=i
100	  if(i.eq.9)go to 110
	    pos=pos+delo
	    call drive(no,pos,0,0,ierr)
c	    if(ith.ne.0)call drive(2,omeg0-(pos-tth0)*0.5,0,0,ierr)
	    call wait
110	continue	  
	if((nstep.le.2).or.(nstep.ge.8))go to 90
C
C CALCULATION OF PEAK BY 9-STEP PARABOLA LEAST-SQUAREs
C
	coefb=sumxy/60.0
	coefa=-sumy/46.2+sumx2y/308.
	if(coefa.ge.0.0)go to 200	!parabola has wrong curvature
	offset=-0.5*coefb/coefa
c next test changed from goto90 by RJA 5.11.93 to prevent infinite  loop
	if(abs(offset).gt.2.5)go to 200 
	posmax=curpos+delo*offset
	call drive(no,posmax,0,0,ierr)
c	if(ith.ne.0)call drive(2,omeg0-(posmax-tth0)*0.5,0,0,ierr)
	call wait
	return
c increment factor for delo changed from 2.0 to 1.5 RJA 5.11.93
c and nloop trap increased from 2 to 3 to compensate
200	delo=delo*1.5
	nloop=nloop+1
	if(nloop.lt.3)go to 90
	rtmax=0			!send fault indicator
	type *,' CENTERING FAILURE IN MTSTEP, MOTOR = ',no,' COEFA = ',coefa
	return
	end
	subroutine refine
	include 'single.inc'
	dimension tphi(3),th(3),angl(4)
c
	if(nlst.gt.0)go to 110
	type 2
2	format(' Reflection List Empty, REFINE Aborted.')
	return
c
c type warning if badcrystal set
c
110	if(badcrystal)type 91
91	format(' FINAL OMEGA SCAN WILL BE FIT BY PARABOLIC METHOD')
c
c start loop over reflections
c
	do 200 i=1,nlst
	  if(iuse(i) .eq. 0)goto200
c
c check for indexing
	  j=abs(ihlst(1,i))+abs(ihlst(2,i))+abs(ihlst(3,i))
	  if(j .eq. 0)then
		do j=1,4
		   cangl(j)=anglst(j,i)
		enddo
		call softlimits(cangl,0,ierr) 
		if(ierr .eq. 1)goto199
	  else
		do j=1,3
	           th(j)=ihlst(j,i)
		enddo
	  	call calphi(th,tphi)
	  	call calca(tphi,cangl,ierr)
	  	if(ierr.eq. 1)goto 199
	  endif
	  call drivep(0,ierr)
	  call wait
	  call center(ifnd,0,width,tthzero)
	  if(ifnd.ne.0)then
   	      write(*,8)i
8	      format(' Reflection No.',i3,' Not Found or too Weak.')
	      goto198
   	  endif
   	  call count(5.0)
   	  call rdallmtr(angl)
	  angl(1)=angl(1)-tthzero
	  call sclwt
	  call ofbeam
	  rate=rdsclr()/rdtime()
	  type 10,i,angl,tthzero,int(rate)
10	  format(' Reflection No.',i3,' at:',4f8.2,' Tthzero ',f5.2,
     1           ' Rate ',i7)
	  do 180 j=1,4
	    anglst(j,i)=angl(j)
180	  continue
	    lstcount(i)=rate
	    iuse(i)=1
	    iref(i)=2
	    call update
	    goto200
199	    type 6,i
6	    format(' Reflection No.',i3,' not Accessible.')
198	    lstcount(i)=0
	    iuse(i)=0
	    iref(i)=0
	    call update
200	continue
	return
	end
	options/extend
	subroutine pklsq(xd,y,ndata,x,inv)
c
c perform fit of Gaussians for alpha1 and alpha2 doublet
c
c xd is array of positions, y is array of counts, ndata is number of points
c x is return array of parameters of equation: BG, GAMMA, TH1, I0
c
c Equation fitted is:
c
c  I(THETA)=I0*DTTH*CG/GAMMA*(EXP(-(AG*(THETA-TH1)/GAMMA)**2)
C	+0.5*EXP(-(AG*(THETA-TH1-CSEP*TAN(2TH/2))/GAMMA)**2))+BG
C
C    WHERE:
C	THETA IS ANGLE OF OBSERVATION
C	I0 IS INTEGRATED INTENSITY OF ALPHA1 PEAK
C	DTTH IS STEP SIZE
C	CG AND AG ARE CONSTANTS TO NORMALIZE INTEGRAL
C	TH1 IS THETA ANGLE OF ALPHA1
C	2TH is two theta for peak.
C	GAMMA IS FULL WIDTH AT HALF MAXIMUM FOR PEAK
C	CSEP IS CONSTANT RELATING SEPARATION OF ALPHA1 AND ALPHA2 
c		(0.345 FOR Mo, 0.141 for Cu)
C
	include 'single.inc'
	dimension xd(1),y(1),x(1)
	dimension a(4,4),bb(4),deriv(4),temp(8),sigx(4),bsol(4)
     1	,sav1(4)
	DATA CG,AG/0.939437279,1.66510922/ 
c
c calculate step size and get tan(two theta/2)
	dtth=abs(xd(2)-xd(1))
	tanth = tan(rdmtr(1)/114.59156)
c
c start cycle - clear matrices
c
	ncyc=0
100	ncyc=ncyc+1
	if(ncyc.gt.100)then
		type 9997
9997		format(' No Convergence in PKLSQ')
		write(98,9997)
		write(98,9999)(x(i),i=1,4)
9999		format(4g14.6)
		write(98,9998)(xd(i),y(i),i=1,ndata)
9998		format(1x,2g14.6)
		inv=1
		return
	endif
	do 110 i=1,4
	  bb(i)=0.0
	  do 110 j=1,4
	    a(j,i)=0.0
110	  continue
120	continue
	gof=0.0
c loop through data
	do 150 n=1,ndata
	  theta=xd(1)+(n-1)*dtth
	  dth1=theta-x(4)
	  dth2=dth1-csep*tanth
	  exp1=exp(-(ag*dth1/x(2))**2)
	  exp2=exp(-(ag*dth2/x(2))**2)
	  se=exp1+0.5*exp2
	  yc=x(3)*dtth*cg*se/x(2)+x(1)
	  delt=y(n)-yc
	  if(y(n) .le. 0.0)then
		weight=1.
	  else
	  	weight=1.0/y(n)
	  endif
c calculate derivatives analytically
	  deriv(1)=1.0
	  deriv(3)=(yc-x(1))/x(3)
	  deriv(2)=x(3)*dtth*cg/x(2)*(-se/x(2)+ag**2/x(2)**3*
     1	    (2.0*dth1**2*exp1+dth2**2*exp2))
	  deriv(4)=x(3)*dtth*cg*ag**2/x(2)**3*(2.0*dth1*exp1
     1	    +dth2*exp2)
c augment least-sqares matrices
	  do 140 i=1,4
	    bb(i)=bb(i)+weight*delt*deriv(i)
	    do 130 j=1,4
	      a(j,i)=a(j,i)+weight*deriv(i)*deriv(j)
130	    continue
140	  continue
	  gof=gof+weight*delt**2
150	continue
	gof=gof/(ndata-4)
c condition matrix
	do i=1,4
c introduce change from LWF 17-nov-1992
	  if(a(i,i) .lt. 1.0e-20)then
	    write(*,'(a)')' Singularity in PKLSQ'
	    inv=1
	    return
	  endif
	  sav1(i)=1.0/sqrt(a(i,i))
	enddo
	do i=1,4
	  do j=1,4
		a(i,j)=a(i,j)*sav1(i)*sav1(j)
	  enddo
	enddo
c invert matrix
	call mxlneq(a,4,4,det,jrank,1.0e-6,bsol,0)
	if(jrank.ne.4)then
   		write(*,'(a)')' Singularity in PKLSQ'
   		inv=1
		return
	endif
c correct for conditioning operations
	do i=1,4
	  do j=1,4
	    a(i,j)=a(i,j)*sav1(i)*sav1(j)
	  enddo
	enddo
c calculate solution vector
	do 170 i=1,4
	  bsol(i)=0.0
	  do 160 j=1,4
	    bsol(i)=bsol(i)+a(i,j)*bb(j)
160	  continue
170	continue
c apply increments to parameters with 0.7 damping factor
	irep=0
	do 190 i=1,4
	  x(i)=x(i)+0.7*bsol(i)
	  if(gof*a(i,i).ge.1.0e-8)go to 180
	    inv=1
	    return
180	  sigx(i)=sqrt(gof*a(i,i))
c check for convergence
	  if(abs(bsol(i)/sigx(i)).gt.0.02)irep=1
190	continue
c	write(98,'(a,i3,4g14.3,a,g12.6)')' Refined Results, Cycle ',ncyc,
c     1  (x(j),j=1,4),' GOF',gof
c loop back for another cycle if not converged
	if(irep.ne.0)go to 100
	inv=0
	return
	end
