C PGM PZSOURCE
	PARAMETER NFFT=1024,NP=NFFT+2,NF=NP/2
	COMPLEX err(nf),zload,zl(3),r(3),c(nf,3),Z0(NF),P0(NF)
	dimension er(4),xlold(3)
	dimension dx(3),xl0(3),xl(3),d(3)
	logical first/.true./,isensw
	data erold/1E9/,dx/3*5E-6/
	data emin/1e-6/,emax/.1/,dmin/.01/,dmax/5./

	CALL INITT(0)
	CALL HEAD("PZSOURCE")
	call sread(nfft,c(1,1),iclk,"c0")
	call sread(nfft,c(1,2),iclk,"C1")
	call sread(nfft,c(1,3),iclk,"c2")

	dt=1e-6*iclk
	fmax=.5/dt
	fmin=fmax/(nf-1)
	fhi=20000.
	flo=500
	accept "flo,fhi=",flo,fhi
	iflo=1+flo*(nf-1)/fmax
	ifhi=1+fhi*(nf-1)/fmax

c set up initital estimates
	xl(3)=xl0(3)=11.85
	xl(2)=xl0(2)=5.94
	xl(1)=xl0(1)=2.75
	d(1)=d(2)=d(3)=7.5  ;dia of tube
	if(isensw(1))accept "xl (mm)=",xl0
	if(isensw(0))accept "d (mm)=",d(1)
	d(2)=d(3)=d(1)
	DO 3 I=1,3
	D(I)=D(I)/1000.
	XL0(I)=XL0(I)/1000.  ;CONVERT TO METER FROM MM
3	CONTINUE
100	continue

	do 15 igrad=1,4  ;compute grad
	call copy(xl0,xl,6)
	if(igrad.ne.1)xl(igrad-1)=xl0(igrad-1)+dx(igrad-1)

c set up loop over freq
	sum=0

	DO 10 I=iflo,ifhi
	FREQ=FMIN*(I-1)

	do 12 j=1,3
	r(j)=c(i,j)
	zl(j)=zload(freq,xl(j),d(j))
12	continue

	call clsqpz(p0(i),z0(i),error,zl,r)
	sum=sum+error
	err(i)=error

10	continue
	sum=sum/(1.+ifhi-iflo)
	er(igrad)=sum
15	continue

	xlen=0
	do 16 i=1,3
	xl(i)=-(er(i+1)-er(1))/dx(i)
	xlen=xlen+xl(i)*xl(i)
16	continue
	xlen=sqrt(xlen)
	type xl0,xlen,er(1)
c	type "grad,len=",xl,xlen
	if(er(1).gt.erold)go to 200

c save for next pass
	erold=er(1)
	call copy(xl0,xlold,6)  ;save last pass

c update length in direction of gradient
	do 17 i=1,3
	xl0(i)=xl0(i)+dx(i)*xl(i)
17	continue
c	go to 100	;try next value

2	continue
	if(.true.)go to 1
	call gwindo(1,3,1,1)
	call dwindo(flo,fhi,dmin,dmax)
	call logtrn(3)

	if(first)call gridl(10,10,0,"freq","p0 & c0")
	if(first)call pltfast(c,nfft,dt)
	CALL PLTFAST(P0,NFFT,DT)

	call gwindo(1,3,1,2)
	call dwindo(flo,fhi,.1,10.)
	if(first)call gridl(10,10,0,"freq","z0")
	call pltfast(z0,nfft,dt)

1	continue
	CALL GWINDO(1,1,1,1)
	call dwindo(flo,fhi,emin,emax)
	call logtrn(3)
	if(first)call gridl(10,10,0,"freq","err")
	call pltfast(err,nfft,dt)

	first=.false.
	call bell
	go to 100

c recompute over entire freq range, and save on disk
200	continue
	do 50 i=1,nf
	freq=fmin*(i-1)
	do 14 j=1,3
	xl(j)=xlold(j)
	r(j)=c(i,j)
	zl(j)=zload(freq,xl(j),d(j))
14	continue

	call clsqpz(p0(i),z0(i),err,zl,r)
50	continue

	call newpag
	CALL HEAD("PZSOURCE")
	type "len'",xl
	call gwindo(1,2,1,1)
	call dwindo(fmin,fmax,.01,5.)
	call logtrn(3)
	call gridl(10,10,0,"freq","p0")
	call pltfast(p0,nfft,dt)

	call gwindo(1,2,1,2)
	call dwindo(fmin,fmax,.2,5.)
	call gridl(10,10,0,"freq","z0")
	call pltfast(z0,nfft,dt)

	open 1,"psource"
	itype=0
	volt=10
	navg=0
	tcor=0
	write(1)itype,nfft,volt,navg,iclk,tcor
	write(1)p0
	close 1

	open 1,"zsource"
	write(1)itype,nfft,volt,navg,iclk,tcor
	write(1)z0
	close 1
	stop done
	end