      options/extend
      subroutine search
c
c routine to search for diffraction maxima by step scanning 2theta and
c    chi with phi slews
c
c extended rja 10-feb-92 to include omega search in fixed-phi mode
c
      implicit none
      include 'single.inc'
      include 'device.inc'
      real*4 phivec(3),phitmp(3),bg(5),chor,cvert,cdist,tthst,tthen,tthpos
     1,chist,chien,chipos,delchi,sgn,phist,phien,delphi,rdsclr,rdtime,rdmtr
     2,rate,phi,cmax,cmin,tempr,tempt,deltth,temp,thresh,phipos,phiend,pkmax
     3,pkloc,time,pos,width,vlen,tmplen,cosv,pktime,stpsiz
      integer*4 ierr,nfnd,nstart,nstp,imotor,nsig,ntemp,irefind,i
     1,ibgpnt,ns,ifnd,j,nstep,npeaks
      byte resp
      data cdist/180.0/,delphi/2.0/,nsig/8/,delchi/1.0/
c
      call prnt_on
1     type 2
2     format('$Search Routine, Enter Number of Refl. to Find: ')
      read(5,3,end=9999,err=1)nfnd
3     format(i10)
      if(nlst.ne.0)then
114	type 4,nlst
4       format('$Reflection List Contains',i3,' Items.'
     1   ,' Do You Wish to Clear the List? ')
	read(5,6,end=9999,err=114)resp
6       format(a1)
	if(resp.eq.'Y'.or.resp.eq.'y')nlst=0
      endif
      nstart=nlst
118   type 8
8     format('$Enter Horizontal and Vertical Counter Aperature in mm: ')
      read(5,10,end=9999,err=118)chor,cvert
10    format(10f10.0)
11    type 12
12    format('$Enter Initial and Final 2Theta: ')
      read(5,10,end=9999,err=11)tthst,tthen
      tthpos=tthst
      call drive(1,tthpos,0,1,ierr)
	if(ierr .eq. 1)then
		type 2040
2040		format('  Starting 2Theta out of bounds')
		goto11
	endif
13    type 14
14    format('$Enter Initial and Final Chi: ')
      read(5,10,end=9999,err=13)chist,chien
      delchi=sgn(chien-chist)
      call drive(3,chist,0,1,ierr)
	if(ierr .eq. 1)then
		type 15
15		format('  Starting Chi OUT-OF-BOUNDS')
		goto13
	endif
c
c branch on fixed_phi
c note all angle arguments still phist, phiend etc even if fixed_phi=1
c 
	if(fixed_phi .eq. 0)then
		imotor=4
1015  		type 16
16    		format('$Enter Initial and Final Phi: ')
      		read(5,10,end=9999,err=1015)phist,phien
      		delphi=delphi*sgn(phien-phist)
      		call drive(4,phist,0,1,ierr)
	        call drive(2,0.0,0,0,ierr)
		if(ierr .eq. 1)then
			type 17
17			format('  Starting Phi OUT-OF-BOUNDS')
			goto1015
		endif
	else
		imotor=2
1020  		type 1025
1025    	format('$Enter Initial and Final Omega: ')
      		read(5,10,end=9999,err=1020)phist,phien
      		delphi=delphi*sgn(phien-phist)
      		call drive(2,phist,0,1,ierr)
	        call drive(4,0.0,0,0,ierr)
		if(ierr .eq. 1)then
			type 1030
1030			format('  Starting Omega OUT-OF-BOUNDS')
			goto1020
		endif
	endif
c                         
c request sig level for search to recognise peak
c
1049	type 1050,nsig
1050	format('$ Input significance level for peak search,',
     1   ' <RETURN> for  ',i2,'?')
	read(5,'(i5)',err=1049)ntemp
	if(ntemp .gt. 0)nsig=ntemp
c
c ask if ref to be done at end of search
c
1059	type 1060
1060	format('$Input 1 for REF at End of Search: ')
	read(*,1070,err=1059)irefind
1070	format(i5)
        deltth=2.0*57.2958*atan(0.5*chor/cdist)
        deltth=aint(10.0*deltth+0.5)*0.1*sgn(tthen-tthst)
	type *,'Delta TTH:',deltth
c Wait for motors done
	call wait
c Enter 2Theta Loop
200   if(fixed_phi .eq. 0)call drive(2,0.0,0,0,ierr)
      call drive(1,tthpos,0,1,ierr)
      if(ierr .eq. 1)goto460
      delchi=57.2958*2.0*atan(0.5*cvert/(cdist*sin(tthpos/57.2958)))
     1    *sgn(delchi)
      delchi=aint(delchi*10.0+0.5)*0.1
      type *,'Delta Chi:',delchi
      chipos=chist
C Enter Chi Loop
300   call drive(3,chipos,0,1,ierr)
      if(ierr .eq. 1)goto455
      call drive(imotor,phist,0,1,ierr)
      call wait
c
c load initial background count
c
      call onbeam
      call count(1.0)
      call sclwt
      temp=rdsclr()
      thresh=(temp+nsig*sqrt(temp))/rdtime()
      do i=1,5
		bg(i)=thresh
      enddo
      thresh=5.*thresh
      ibgpnt=1
      phipos=phist
c
C Enter Phi Loop
c
400   phiend=phipos+delphi
      call count(100.)
      call drive(imotor,phiend,0,1,ierr)
      if(ierr .eq. 1)goto400
      call wtmtr(imotor)
      call sclhlt
      tempr=rdsclr()
      tempt=rdtime()
      rate=tempr/tempt
      if(rate.lt.thresh/5.)then
c no peak - update rolling average threshold
	  thresh=thresh-bg(ibgpnt)
          bg(ibgpnt)=(tempr+nsig*sqrt(tempr))/tempt
	  thresh=thresh+bg(ibgpnt)
	  ibgpnt=mod(ibgpnt,5)+1
	  go to 450
      endif
C Peak Found - Find Maximum on slewing drive
      pkmax=0.0
      nstp=abs(4.0*delphi)
      do i=1,nstp+1
	phi=phiend-(i-1)*delphi/nstp
	call drive(imotor,phi,0,1,ierr)
	call wait
	call count(0.25)
	call sclwt
	rate=rdsclr()/rdtime()
	if(rate.gt.pkmax)then
	  pkmax=rate
	  pkloc=rdmtr(imotor)
	endif
      enddo
c
c drive back to maximum
c
      call drive(imotor,pkloc,0,1,ierr)
      call rdallmtr(cangl)
      if(pkmax.lt.30.0)then
c
c here if peak too weak to centre
c
	type 401,cangl,int(pkmax)
401	format(7x,'Weak Maximum Found at     ',4f8.1,'Rate: ',i2,' cps')
      else
c
c here if peak strong enough to center 		
c
c -initial omega and check for peak again
	call mtstep(2,.10,cmax,cmin,.5,ns,ierr)
	if(ierr .eq. 1)goto430
	if(cmax .lt. minct)goto430
	if(cmin/cmax .gt. 1./maxmin)goto430
c
c theta-2theta and check for peak
	call mtstep(1,.2,cmax,cmin,.5,ns,ierr)
	if(cmax .lt. minct)goto430
	if(cmin/cmax .gt. 1./maxmin)goto430
	if(ierr .eq. 1)goto430
c
c omega again, then chi
	call mtstep(2,.05,cmax,cmin,.5,ns,ierr)
	if(ierr .eq. 1)goto430
	call mtstep(3,.3,cmax,cmin,.5,ns,ierr)
	if(ierr .eq. 1)goto430
	if(cmax .lt. minct)goto430
c
c do and fit omega scan
c
	time=555.5/cmax
	if(time .lt. 0.25)time=0.25
	if(time .gt. tmax)time=tmax
	call omstep(ifnd,pos,time,width)
	if(ifnd .ne. 0)goto430
	call drive(2,pos,0,1,ierr)
	if(ierr .eq. 1)goto430
c
c success at centering-check against list
c
        call rdallmtr(cangl)
	if(nlst.ne.0)then
	  call calch(cangl,phivec)
	  vlen=phivec(1)**2+phivec(2)**2+phivec(3)**2
	  do i=1,nlst
	    call calch(anglst(1,i),phitmp)
	    tmplen=phitmp(1)**2+phitmp(2)**2+phitmp(3)**2
	    if(abs(vlen-tmplen).le.0.005)then
c vector length matches one of the list items - check vector components
	      cosv=(phitmp(1)*phivec(1)+phitmp(2)*phivec(2)
     1        +phitmp(3)*phivec(3))/sqrt(vlen*tmplen)
	      if(abs(cosv).gt.0.999)then	
c vectors are parallel
		type 402,i,cangl
402		format(7x,'Duplicate  No.',i3,', Angles:',4f8.2)
	        goto440
	      endif
	    endif
	  enddo
	endif
c Vector is Unique: LOAD INTO LIST
	nlst=nlst+1
c clear indices
	do j=1,3
	  ihlst(j,i) = 0.0
	enddo
	do j=1,4
	  anglst(j,nlst)=cangl(j)
	enddo
	iref(nlst)=2
	iuse(nlst)=1
	call onbeam
	call count(1.0)
	call sclwt
	rate=rdsclr()/rdtime()
	lstcount(nlst)=rate
	type 18,nlst,cangl,int(rate),width
18      format(/7x,'New Peak,  No.',i3,', Angles:',4f8.2/
     1  7x,'Rate:',i8,' cps,  Width:',8x,f8.2)
	pktime=555.5/rate
	if(pktime .lt. 0.5)pktime=0.5
	if(pktime .gt. 3.0)pktime=3.0
	stpsiz=0.025
c
c set number of steps to 5*full width, subject to total scan
c being less than 1.5 deg wide
c
	width=2.*width	
	nstep=5.0*width/stpsiz+0.5
	if(nstep .gt. 60)nstep=60
	call pkprof(nstep,stpsiz,pktime,2,ierr)     !generate peak profile on printer
	call update
	if(nlst-nstart.ge.nfnd)go to 500
      endif
      goto440
c
c failure of peak centering
430   type 431,cangl,int(pkmax)
431   format(7x,'Peak Centering Failed at:',4f8.1,' Rate:',i6,' cps')
c
c go to end of next step and restart search
440   call drive(1,tthpos,0,0,ierr)
      if(fixed_phi .eq. 0)call drive(2,0.0,0,0,ierr)
      call drive(3,chipos,0,0,ierr)
      phiend=phiend+delphi
      call drive(imotor,phiend,0,0,ierr)
      call wait
      call onbeam
c 
c increment phi 
c
450   phipos=phiend
      if((phiend-phist)*(phien-phiend).gt.0.0)go to 400  !loop for more Phi's
c End of Phi Scan - Advance to Next Chi
455   delphi=-delphi
      temp=phist
      phist=phien
      phien=temp
      chipos=chipos+delchi
      if((chipos-chist)*(chien-chipos).gt.0.0)go to 300
C End of Chi Steps - Advance to next 2Theta
460   delchi=-delchi
      temp=chist
      chist=chien
      chien=temp
      tthpos=tthpos+deltth
      if((tthpos-tthst)*(tthen-tthpos).gt.0.0)go to 200
C End of 2Theta Steps - Done with Problem
500   call ofbeam
c
c drive to zero
      do i=1,4
		cangl(i)=0.
      enddo
      call drivep(0,ierr)
      npeaks=nlst-nstart
      type 30,npeaks
30    format(/' End of Search:',i3,' Peaks Found.')
      if(irefind .eq. 1)call refine
      call update
c     call prnt_off
9999  return
      end
c
c
c
	subroutine cone
c
c does cone search around diffraction vector
c
	implicit none
	include 'single.inc'
	real*4 pmat(3,3),rdelt(3,3),rgamma(3,3),angl(4),ttht
     1	,rprime(3,3),rtemp(3,3),rad,hphit(3),bg(5),tth,q,cnangl
     2  ,rdmtr,rdsclr,rdtime,gamma,temp,tempr,tempt,rate,thresh
	integer*4 nsig,ifirst,i,ierr,ibgpnt
	character*1 iresp
	data rad/57.29578/
c
c Print initial info
c
	call wait
5	type 10
10	format(/'$Cone Search - Is Diffractometer at a Reflection?')
	read(5,'(a)',end=9999,err=5)iresp
	if(iresp.ne.'Y'.and.iresp.ne.'y')then
	    type 20
20	    format(' Position Diffractometer and Retry CONE')
	    return
	endif
c
c Get current angles
c
	call rdallmtr(angl)
c Calculate pmat
	call psimat(angl,pmat)
c Get 2theta and cone angle
1000	type 30
30	format(/'$Enter Cone Angle to Second Reflection and Its 2Theta? ')
	read(5,'(2f10.0)',end=9999,err=1000)cnangl,tth
	if(tth.le.0.0)tth=angl(1)
	q=2.0*sin(0.5*tth/rad)/wavel
	if(cnangl.le.0.0)then
c Impossible cone angle
	    type 40
40	    format(' CONE Angle too small - Try again.')
	    go to 1000
	endif
c Calculate R(DEL)T
	rdelt(1,1)=cos(cnangl/rad)
	rdelt(1,2)=0.0
	rdelt(1,3)=-sin(cnangl/rad)
	rdelt(2,1)=0.0
	rdelt(2,2)=1.0
	rdelt(2,3)=0.0
	rdelt(3,1)=sin(cnangl/rad)
	rdelt(3,2)=0.0
	rdelt(3,3)=cos(cnangl/rad)
c Get Significance Number
49	type 50
50	format(/'$Type Significance Number <DEF=8>?')
	read(5,'(i10)',err=49)nsig
	if(nsig.le.0)nsig=8
	type 52
52	format(/' Locations of Significant Intensity'/)
c
c Start Loop for Rotation about the Cone
c
	ifirst=0
	do gamma = 0.,358.,2.
c Form R(GAMMA)T
	    rgamma(1,1)=1.0
	    rgamma(2,1)=0.0
	    rgamma(3,1)=0.0
	    rgamma(1,2)=0.0
	    rgamma(2,2)=cos(gamma/rad)
	    rgamma(3,2)=sin(gamma/rad)
	    rgamma(1,3)=0.0
	    rgamma(2,3)=-sin(gamma/rad)
	    rgamma(3,3)=cos(gamma/rad)
c Calculate R'
	    call matmul(rdelt,rgamma,rtemp)
	    call matmul(rtemp,pmat,rprime)
c Get h(phi) for New orientation
	    do i=1,3
		hphit(i)=rprime(1,i)*q
	    enddo
c Get Angles
	    call calca(hphit,cangl,ierr)
	    if(ierr.eq.0)then
c Angles OK - Try to drive to them
		call drivep(0,ierr)
		if(ierr.eq.0)then
		    call wait
		    if(ifirst.eq.0)then
			ifirst=1
c
c load initial background count
c
			call onbeam
			call count(1.0)
			call sclwt
			temp=rdsclr()
			thresh=(temp+nsig*sqrt(temp))/rdtime()
			do i=1,5
			    bg(i)=thresh
			enddo
			thresh=5.*thresh
			ibgpnt=1
		    endif
c Do a Count
		    call count(0.5)
		    call sclwt
		    tempr=rdsclr()
		    tempt=rdtime()
		    rate=tempr/tempt
		    if(rate.lt.thresh/5.)then
c no peak - update rolling average threshold
			thresh=thresh-bg(ibgpnt)
			bg(ibgpnt)=(tempr+nsig*sqrt(tempr))/tempt
			thresh=thresh+bg(ibgpnt)
			ibgpnt=mod(ibgpnt,5)+1
		    else
			type 60,cangl,int(rate)
60			format(' Angles: ',4f8.2,' Rate: ',i8)
		    endif
                endif
	    endif
	enddo
	call ofbeam
	do i=1,4
	    cangl(i)=0.0
	enddo
	call drivep(0,ierr)
9999	return
	end
c
c
	SUBROUTINE ANGLETABLE
c
c subroutine to produce table of inter-reflection angles
c
	implicit none
	include 'single.inc'
	real vangle(13),phivec(3),phitmp(3),vlen,tmplen,cosv
	character*255 fstring
	integer*4 num,i,n_list,ipt,j,k,istart
	integer*4 listpnt(listlen)
	byte esc/27/
c
c Fill in pointer list
c
	n_list=0
	do i=1,nlst
	  if(iuse(i) .ne. 0)then
	    n_list=n_list+1
	    listpnt(n_list)=i
	  endif
	enddo
c
c Start loop through list
c
	do i=1,n_list,12
	  num = min(12,n_list-i)
c
c Write Header
c
	  write(*,'(/4x,12i6)')(listpnt(j),j=i,i+num-1)
c
c loop over reflections in this strip
c
  	  do j=i+1,n_list
	    fstring='(1x,i3,1h:,'
	    istart=11
	    call calch(anglst(1,listpnt(j)),phivec)
	    vlen=phivec(1)**2+phivec(2)**2+phivec(3)**2
	    ipt=0
	    do k=i,min(j-1,12+i-1)
	      call calch(anglst(1,listpnt(k)),phitmp)
	      ipt=ipt+1
	      tmplen=phitmp(1)**2+phitmp(2)**2+phitmp(3)**2
c
c calculate angle
c	      
	      cosv=(phitmp(1)*phivec(1)+phitmp(2)*phivec(2)
     1        +phitmp(3)*phivec(3))/sqrt(vlen*tmplen)
	      vangle(ipt)=57.2958*acos(cosv)
	      if(vangle(ipt) .lt. 10.)then
		fstring=fstring(1:istart)//'3x,a1,3h[7m,f3.1,'
	        istart=istart+17
	      else
	 	fstring=fstring(1:istart)//'a1,3h[0m,f6.1,'
		istart=istart+14
	      endif
	    enddo
c
c end of one line
c
	    fstring=fstring(1:istart)//'a1,3h[0m)'
	    write(*,fstring)listpnt(j),(esc,vangle(k),k=1,ipt),esc
c
	  enddo
	enddo
	return
	end
