c HL modifications (Nov. 26, 1996)
c    (1) w(ii) set to 1. everywhere w(ii) initially .ge. 10., not
c        just near edge
c    (2) cor2(i) set to zero unless wav and cnorm both .gt. 0.
c    (3) below line 336: changed cx(i,itpeak) to cx(i,ihpeak)
c HL modifications (Apr. 17, 1997)
c    set ztrust = 0.6
c HL modifications (May 16, 1997)
c    declare errmsg as character*80
c HL modifications (Nov 7, 1997)
c    hardwire w0g = 3.643 and dw = 3.e-4 because of problems
c    with reading parameters at end of long iraf headers
c

	program	corel
	dimension gspec(2000),tspec(2000)
	dimension w(2000)
	real*8 ww,wav
	dimension cor(4000),ip(10),rr(20),pmax(20)
	dimension tnorm(4000),cx(4000,20),cor2(4000)
	dimension fcor(4000),fcor2(4000),iflag(20)
	integer	axleng(7),axlent(7)
        character*80 errmsg
c
c	Get images names
c	G means	the data spectrum and T	means the template spectrum
c
	inqpeak=1
	rr(0)=0.0
	ip(0)=0.0
	call imopen('finr',1,img,ier)
	if(ier.ne.0) goto 91
	write(*,*) 'open galaxy'
	call imopen('template',1,imt,ier)
	if(ier.ne.0) goto 91
	write(*,*) 'open template'
c
	call imgsiz(img,axleng,naxis,dtype,ier)
	if(ier.ne.0) goto 91
	write(*,*) axleng(1), axleng(2) 
	call imgsiz(imt,axlent,naxist,dtype,ier)
	if(ier.ne.0) goto 91
	write(*,*) axlent(1), axlent(2) 
c	axleng and axlent are the axis lengths of the spectrum and temp
cc	call imgkwr(img,"w0",w0g,ier)
cc	if(ier.ne.0) goto 91
        w0g = 3.643
	call imgkwr(imt,"w0",w0t,ier)
	write(*,*) w0g,	w0t
c	w0t and	w0g are	the wavelength zer0points of the two spectra
	if(ier.ne.0) goto 91
cc	call imgkwr(img,"wpc",dw,ier)
cc	if(ier.ne.0) goto 91
        dw = 3.e-4
	call imgkwr(imt,"wpc",dw2,ier)
	if(ier.ne.0) goto 91
	if(dw.ne.dw2) goto 90
	write(*,*) dw, dw2
c	this program will not work unless the two spectra have the same
c	(log lambda)/pixel
c	call imgkwr(imt,"v0",v0,ier)
c	if(ier.eq.0) goto 92 
	v0=0.
c	v0=cz for the galaxy template 
c	for multimake, corrected afterwards
c92	write(*,*) 'v0 (template) ',v0, v0/2.9989e5
	naxleng=axleng(1)
	naxlent=axlent(1)
	nspectra=axleng(2)
	ntemp=axlent(2)
	dlambda=(w0g-w0t)/dw
	ndlam=int(dlambda)
	errlam=(dlambda-float(ndlam))*dw
	write(*,*) dlambda, ndlam, errlam
c	naxleng and naxlent are the lengths of	the two	spectra
c	dlambda	is the difference in starting wavelengths
c	ndlam is the difference	in integer pixels
c	errlam is the remainder	is the difference
c	
	narray=naxleng+naxlent-1
c
	izp=naxlent-ndlam
	write(*,*) 'izp= ',izp
	minz=izp
	naxlenc=int(0.3/dw)
	maxz=naxlenc+izp
c	max z=1.0
c	outputspectrum is from z=0 to z=1.0
c	uncomment below to output whole correlation
c	naxlenc=narray
c	minz=1
c	maxz=narray  
c
	axleng(1)=naxlenc
	axleng(2)=axleng(2)
	if(ntemp.gt.1) naxis=3
	if(naxis.eq.3) axleng(3)=ntemp 
	write(*,*) naxis,axleng(1), axleng(2), axleng(3)
	call imcrea("cor",axleng,naxis,dtype,ier)
	if(ier.ne.o) goto 91
	write(*,*) 'create cor'
	call imopen("cor",3,imc,ier)
	if(ier.ne.o) goto 91
	call imakwi(imc,"dispaxis",1,"",ier)
	if(ier.ne.0) goto 91
	write(*,*) 'open cor'
c
c	ztrust=(0.5*naxleng*dw+w0g)-w0t
c	fudge below!	
ccc	ztrust=(0.5*naxleng*dw+w0g)-3.56
ccc	ztrust=10**(ztrust) - 1
        ztrust = 0.6
	write(*,*) 'ztrust= ', ztrust
	rlimit=min(500,0.15/dw)
	itrust=izp+alog10(ztrust+1)/dw
	write(*,*) itrust, narray
	rlimit=min(rlimit,narray-itrust)
	write(*,*) 'rlimit= ', rlimit 

c
	do 9999 iap=1,nspectra 
	code=0

	call imgl3r(img,gspec,iap,1,ier)
	if(ier.ne.0) goto 91
	call imgl3r(img,w,iap,2,ier)
	if(ier.ne.0) goto 91
c
	ns=0
	tsig=0
	gsig=0
	wscale=0
	do 600 ii=1,naxleng
	gsig=gsig+(gspec(ii)*gspec(ii))
	wscale=wscale+w(ii)
600	continue
	wscale=wscale/float(naxleng)
c	write(*,*) 'wscale= ', wscale
	do 603 ii=1,naxleng
	w(ii)=w(ii)/wscale
	edge=min(ii,naxleng-ii)
	if(w(ii).le.0) then
		w(ii)=0.
	else  
		if(edge.le.10) then
			w(ii)=0.2
		else
			w(ii)=1./(w(ii)**2)
		endif
	endif
cc	if(w(ii).ge.10..and.edge.le.40.)  w(ii)=1.0
        if (w(ii) .ge. 10.) w(ii)=1.0
603	continue
c
	do 1100 intemp=1,ntemp	
	call imgl2r(imt,tspec,intemp,ier)
	if(ier.ne.0) goto 91
c
	do 700 ii=1,naxlent
700	tsig=tsig+(tspec(ii)*tspec(ii))
	cnorm=sqrt(tsig*gsig)
c
	tspec(0)=0.0
	gspec(0)=0.0
c
c	i counts the pixel of the correlation spectrum
	do 100 i=1,narray   
	cor(i)=0.0
	cor2(i)=0.0
	tnorm(i)=0.0
	npts=0      
	tn=0.
	gn=0.
	wsum=0.
c	n counts the length of the data	spectrum
c	npts counts the number of points actually used
c
	do 200 n=1,naxleng
	kdl=naxlent-i+n
c	x-correlation is shift,	multiply and add
c	kdl is the shift between the spectra in pixel units
	if(kdl.gt.naxlent) goto 200
	if(kdl.le.0) goto 200
	if(tspec(kdl).eq.0.) goto 200
	if(gspec(n).eq.0.) goto 200
c	don't include in the sum any points where
c	you run	off the	ends of the template spectrum
c	points where either one	of the spectra are equal to zero
c	relative weights are 1/sigma**2
	if(w(n).eq.0) goto 200 
	ww=w(n)
	cor(i)= cor(i) + gspec(n)*tspec(kdl)*ww
	npts=npts+1
	gn= gn + gspec(n)*gspec(n)*ww
	tn= tn + tspec(kdl)*tspec(kdl)   
	wsum=wsum+ww
200	continue
	if(npts.eq.0) goto 100
	wav=wsum/npts
	cor(i)=(cor(i)/cnorm)/wav
	tnorm(i)=sqrt(tn*(gn/wav))/cnorm
100	continue
c
	npeaks=0
	do 998 ipix=minz,maxz  
	if(cor(ipix).lt.cor(ipix-1).or.cor(ipix).le.cor(ipix+1)) goto 998 
	if(cor(ipix-1).eq.0.or.cor(ipix+1).eq.0) goto 998
	imax=ipix
	x0=(cor(imax+1)-cor(imax-1))
	x0=x0/(2*(2*cor(imax)-cor(imax-1)- cor(imax+1)))
c
c	get hmax
	hmax=cor(imax)
	hn=hmax/tnorm(imax)
c
	call rxcor(imax,cor,narray,iedge,rlimit,intemp,hmax,r)
	if(r.lt.2.0) goto 998
c
	npeaks=npeaks+1
	ip(npeaks)=imax
c
c	recalculate cor removing aliassing
        iflag(npeaks)=0
        if(iedge.eq.1.and.intemp.eq.5) iflag(npeaks)=1
c
	do 601 i=1,narray   
	cor2(i)=0.0
	wsum=0.
	cx(i,npeaks)=0.0
	id=imax-i 
	npts2=0
	do 602 n=1,naxleng
	kdl=naxlent-i+n
	if(kdl.le.0) goto 602
	if(kdl.gt.naxlent) goto 602
	if(gspec(n).eq.0) goto 602
	if(tspec(kdl).eq.0) goto 602
c
	if((n+id).le.0) goto 602
	if((n+id).gt.naxleng) goto 602
	if(gspec(n+id).eq.0) goto 602
	if((n+ndlam+izp-imax).le.0) goto 602
	if((n+ndlam+izp-imax).gt.naxlent) goto 602
	if(tspec(n+ndlam+izp-imax).eq.0) goto 602
	if(w(n).le.0) goto 602 
	ww=w(n)
	wsum=wsum+ww
	cor2(i)= cor2(i) + gspec(n)*tspec(kdl)*ww
	npts2=npts2+1
602	continue
	wav=wsum/npts2

	if (wav .gt. 0. .and. cnorm .gt. 0.) then
           cor2(i)=cor2(i)/wav/cnorm
        else 
           cor2(i) = 0.
        end if

cc	cor2(i)=cor2(i)/wav	
cc	cor2(i)=cor2(i)/cnorm
        if(iflag(npeaks).eq.1) cor2(i)=cor(i) 
	cx(i,npeaks)=cor2(i)
601	continue
c
	x0=(cor2(imax+1)-cor2(imax-1))
	x0=x0/(2*(2*cor2(imax)-cor2(imax-1)- cor2(imax+1)))
c
	call rxcor(imax,cor2,narray,iedge,rlimit,intemp,hmax,r)
	pmax(npeaks)=hmax
	rr(npeaks)=r
c
c	write(*,102) iap,npeaks,imax, x0, hmax,hn,sig,r
102	format(1x,i4,i3,1x,i4,1x,f7.4,1x,f7.4,1x,f7.4,1x,f7.4,1x,f8.2)
c
998	continue
c
 	if (npeaks.eq.0) goto 348
c
355	rmax=0
	pmaxn=0
	hmaxn=0
	do 331 ipk=1,npeaks
	if(rr(ipk).gt.rmax) itpeak=ipk 
	if(rr(ipk).gt.rmax) rmax=rr(ipk)
	hm=pmax(ipk)/tnorm(ip(ipk))
	pm=pmax(ipk)
	if(pm.gt.pmaxn) ihpeak=ipk 
	if(pm.gt.pmaxn) pmaxn=pm
	if(hm.gt.hmaxn) ippeak=ipk 
	if(hm.gt.hmaxn) hmaxn=pm
331	continue  
333	itrust=alog10(1+ztrust)/dw + izp
	rd=(rr(itpeak)*pmax(ihpeak))/(rr(ihpeak)*pmax(itpeak))
	if(ihpeak.eq.itpeak.or.ihpeak.eq.ippeak) goto 336
	if(tnorm(ip(ippeak))/tnorm(ip(ihpeak)).ge.0.8) ihpeak=ippeak  
	if(ihpeak.eq.itpeak.or.ihpeak.eq.ippeak) goto 336
	if(rd.ge.1.0) ihpeak=itpeak 
	if(ip(ihpeak).le.itrust.or.iflag(ihpeak).eq.1)goto 336
	rr(ihpeak)=1.0
	goto 355
c
336	do 345 i=1, narray
	cor(i)=cor(i)/tnorm(ip(ihpeak))
	cor2(i)=cx(i,ihpeak)/tnorm(ip(ihpeak))
cc	cor2(i)=cx(i,itpeak)/tnorm(ip(ihpeak))
345	continue
c
	imax=ip(ihpeak)
c	write(*,*) 'Chose Peak: ', ihpeak
c
	x0=(cor2(imax+1)-cor2(imax-1))
	x0=x0/(2*(2*cor2(imax)-cor2(imax-1)- cor2(imax+1)))
	sx0=sqrt(x0*x0+1)/rr(ihpeak)
	x0=(x0+imax-izp)
	w0=errlam+alog10(1+v0/2.9989e5)
	z=10**(dw*x0+w0)-1
	dz=10**(dw*(x0+sx0)+w0)-z-1
	dv=dz*3e5
	dv=(sqrt(dv*dv + 65.0*65.0))/(1+z)
	w0=w0+(minz-izp)*dw
	write(*,101) iap, x0, sx0, z, dz, dv,  rr(ihpeak), code, intemp
	write(2,101) iap, x0, sx0, z, dz, dv,  rr(ihpeak), code, intemp
101	format(1x,i4,2f7.2,1x,2f8.5,3x,f7.0,3x,f8.2,f5.1,i5)
c
c	kludge OIII redshifts	
	if(intemp.eq.5.and.z.gt.0.34) then
		x0=x0-int(alog10(5007./3727.)/dw)+izp
		imax=int(x0)
		hmax=0
		ixmax=izp
		do 512 io3=imax-8,imax+8
		if(cor(io3).gt.hmax) ixmax=io3
		if(cor(io3).gt.hmax) hmax=cor(io3)
512	continue
		imax=ixmax
		x0=(cor(imax+1)-cor(imax-1))
		x0=x0/(2*(2*cor(imax)-cor(imax-1)-cor(imax+1)))
c
		r=rr(ihpeak)
		sx0=sqrt(x0*x0+1)/r
        	x0=(x0+imax-izp)
        	z=10**(dw*x0+w0)-1
        	dz=10**(dw*(x0+sx0)+w0)-z-1
        	dv=dz*3e5
        	dv=(sqrt(dv*dv + 65.0*65.0))/(1+z)
		write(3,101) iap, x0, sx0, z, dz, dv,  r, 0, 5
	elseif(intemp.eq.5) then 
		write(3,101) iap, 0.,0.,0.,0.,0.,0.,-1., 5
	endif

	goto 346
c	no peak chosen
348	write(2,101) iap, 0.,0.0,0.0,0.0,0.0,0.0,-1.0,intemp
	if(intemp.eq.5) then 
		write(3,101) iap, 0.,0.,0.,0.,0.,0.,-1., 5
	endif
	write(*,*) iap
c
c
346	do 444 if=1,naxlenc
	fcor(if)=cor(if+minz-1)	
	fcor2(if)=cor2(if+minz-1)	
444	continue
c
	call impl3r(imc,fcor,iap,intemp,ier)
	if(ier.ne.0) goto 91
c

1100	continue
	write(*,*)
c
9999	continue
c
	call imakwi(imc,"crpix1",1,"",ier)
	if(ier.ne.0) goto 91
	call imakwr(imc,"w0",w0," ",ier)
	if(ier.ne.0) goto 91
	call imakwr(imc,"wpc",dw," ",ier)
	if(ier.ne.0) goto 91
	call imakwr(imc,"crval1",w0," ",ier)
	if(ier.ne.0) goto 91
	call imakwr(imc,"cdelt1",dw,"",ier)
	if(ier.ne.0) goto 91
	call imakwc(imc,"ctype1",'LINEAR',"",ier)
	if(ier.ne.0) goto 91
	write(*,*) 'write new files'
c
	call imclos(imt,ier)
	if(ier.ne.0) goto 91
	call imclos(img,ier)
	if(ier.ne.0) goto 91
	call imclos(imc,ier)
	if(ier.ne.0) goto 91
c	call imclos(imc2,ier)
	if(ier.ne.0) goto 91
c
	
	goto 777
91	call imemsg(ier,errmsg)
	write(*, '('' Error: '', a80)') errmsg
	goto 777
90	write(*,*) 'Error: files have different dispersion!'
777	continue
	stop
	end
	
	subroutine rxcor(imax,cor,narray,iedge,rlimit,intemp,hmax,r)
	real cor(4000)
c	Calculate sigma of asymmetric part
	r=0.
	sig=0.
	nsum=min(imax,narray-imax-1) 
        iedge=0
        if(nsum.lt.rlimit) iedge=1
        if(nsum.lt.rlimit.and.intemp.ne.5) goto 598
	nsum=rlimit
c 	calculate average
	asav=0
	ns=0
	nnsum=-1*nsum
	do 801 isig=nnsum,nsum
	if(cor(imax+isig).eq.0.or.cor(imax-isig).eq.0) goto 801
	asav=asav+(cor(imax+isig)-cor(imax-isig))
	ns=ns+1
801	continue
	if(ns.lt.50) goto 598
	asav=asav/ns
c	calculate standard deviation
	ns=0
	do 800 isig=nnsum,nsum
	if(cor(imax+isig).eq.0.or.cor(imax-isig).eq.0) goto 800
	sig=sig+(cor(imax+isig)-cor(imax-isig)-asav)**2
	ns=ns+1
800	continue
	sig=sig/(ns-1)
	sig=sqrt(sig)
c	if(hmax.lt.0.000001.and.sig.lt.0.000001) goto 598
	r=hmax/sig
598	return
	end