c
c   searches a MAR IP image for the (110) ring of an orthorhombic
c   form of annealed polyethylene.
c   it then crudely refines the center and xtal-to-plate distance.
c
c                                               Hudel 2/12/93
c
	integer*2 ibuf(2000,2000)
c	integer*4 i4buf(2000,2000)
	real*4 lambda, integ
	character*128 filename, defname
	common ibuf
c
	radius=1000
	diameter=2*radius
	defname = '*.image'
	dtr = 3.14159265358 / 180.0
 	write(6, 10)
10	format('$Image file: ')
	read(5, 20) filename
20	format(a)
c
	init_size = radius
	size_i = radius
	open(
     *		unit=1, status='unknown',
     *		name=filename, defaultfile=defname,
     *		initialsize=init_size, recl=size_i,
     *		extendsize=(size_i*4+511)/512,
     *		access='direct', associatevariable=k_rec,
     *		readonly)
c
	do i1 = 1, diameter
	   read(1, rec=i1+1) (ibuf(i1, i2), i2=1,diameter)
c 	   do i2 = 1, diameter
c               ibuf(i1, i2) = iexpand(ibuf(i1, i2))/2
c 	   enddo
	write(6,*)i1
	enddo
c
1000	continue
	write(*, 24)
24	format('$Wavelength [A]: ')
	read(*, *) lambda
c	lambda = 1.5417
c	lambda = 1.08
 	write(6, 30)
30	format('$Approximate xtal-to-plate distance in mm: ')
	read (*, *) dist

1111	write(*, 25)
25	format('$Resoln of ring: ')
	read(*, *,end=9999,err=9999) res
	if (res.eq.0) goto 9999
c   guessed index of strongest ring for orthorhombic form of
c   IBM (Tom Russel) polyethylene
C	h = 1
C	k = 1
C	l = 0
c   unit cell dims for orthorhombic form of IBM (Tom Russel) polyethylene
C	a = 7.4
C	b = 4.93
C	c = 2.534
c lead nitrate (111) circle
c PbNO3	n = 3
c	res = 1.0 / sqrt((h/a)**2 + (k/b)**2 + (l/c)**2)
c PbNO3	res = 7.8404 / sqrt(real(n))
	type *, 'Searching for circle at ', res, 'A'
	xmm = 2.0 * dist * tan (asin(lambda/res) / 2)
	type *, 'Searching near ', xmm, 'mm from center...'
	xpix = xmm / 0.150
	type *, 'Searching near ', xpix, ' pixels from center...'
c
c	read(*, *) itest
c	if (itest .ne. 0) goto 1111
c
	isumx = 0
	icntx = 0
	isumx = 0
	icntx = 0
	do i = 1, diameter
	   isumx = isumx + ibuf(i, radius)
	   if (ibuf(i, radius) .gt. 0) icntx = icntx + 1
	   isumy = isumy + ibuf(radius, i)
	   if (ibuf(radius, i) .gt. 0) icnty = icnty + 1
	enddo
	avgx = float(isumx) / icntx
	avgy = float(isumy) / icnty
	type *, 'Average intensity along x, y: ', avgx, avgy
c
	maxl = 0
	maxu = 0
	maxr = 0
	maxd = 0
 	do i = xpix - 10, xpix + 10
	   if (ibuf(radius-i, radius) .gt. maxl) then
	      maxl = ibuf(radius-i, radius)
	      maxlx = radius - i
	   endif
	   if (ibuf(radius, radius+i) .gt. maxu) then
	      maxu = ibuf(radius, radius+i)
	      maxuy = radius + i
	   endif
	   if (ibuf(radius+i, radius) .gt. maxr) then
	      maxr = ibuf(radius+i, radius)
	      maxrx = radius + i
	   endif
	   if (ibuf(radius, radius-i) .gt. maxd) then
	      maxd = ibuf(radius, radius-i)
	      maxdy = radius - i
	   endif
 	enddo
	if (maxl .lt. 3*avgx) type *, '*** weak or no ring on left ***',maxl
	if (maxr .lt. 3*avgx) type *, '*** weak or no ring on right ***',maxr
	if (maxu .lt. 3*avgy) type *, '*** weak or no ring up ***',maxu
	if (maxd .lt. 3*avgy) type *, '*** weak or no ring down ***',maxd
	sum1 = 0.0
	sum2 = 0.0
	do i = maxlx-3, maxlx+3
	do j = -5,5
	   sum1 = sum1  +  i * ibuf(i, radius+j)
	   sum2 = sum2  +  ibuf(i, radius+j)
	enddo
	enddo
	rmaxl = sum1 / sum2
	sum1 = 0.0
	sum2 = 0.0
	do i = maxrx-3, maxrx+3
	do j = -5,5
	   sum1 = sum1  +  i * ibuf(i, radius+j)
	   sum2 = sum2  +  ibuf(i, radius+j)
	enddo
	enddo
	rmaxr = sum1 / sum2
	sum1 = 0.0
	sum2 = 0.0
	do i = maxuy-3, maxuy+3
	do j = -5,5
	   sum1 = sum1  +  i * ibuf(radius+j, i)
	   sum2 = sum2  +  ibuf(radius+j, i)
	enddo
	enddo
	rmaxu = sum1 / sum2
	sum1 = 0.0
	sum2 = 0.0
	do i = maxdy-3, maxdy+3
	do j = -5,5
	   sum1 = sum1  +  i * ibuf(radius+j, i)
	   sum2 = sum2  +  ibuf(radius+j, i)
	enddo
	enddo
	rmaxd = sum1 / sum2
c	rxx = (maxrx + maxlx) / 2.0
	rxx = (rmaxr + rmaxl) / 2.0
	ixx = (rxx + 0.5)
c
	ryy = (rmaxu + rmaxd) / 2.0
	iyy = (ryy + 0.5)
	type *, 'Refined center [pixels] at ', rxx, ryy
	type *, 'Refined center [mm] at ', 0.15*rxx, 0.15*ryy
c
c
	diax = rmaxr - rmaxl
	diay = rmaxu - rmaxd
	rrad = (diax + diay) / 4.0
	if (abs(diax-diay) .gt. 2) then	
	   type *, 'Diameters differ by more than 2 pixels!'
	   type *, '   diax, diay: ', diax, diay
	endif
	xdist = 0.150 * (diax + diay) / (8 * tan(asin(lambda/res) / 2))
c	xdist = 0.150 * (diax + diay) / (8 * tan(asin(lambda/(2*res)) )  !is what I think it should be.
	type *, 'Refined distance [mm]: ', xdist
	xres = lambda/(2*sin(.5*atan(0.15*rrad/dist)))
	type *, 'Refined resolution (dist=',dist,')=',xres
c	resmax = lambda / sin(2*atan(0.5*xmm/dist))
c
	vmin = 1e6
	vmax = -1e6
	do ian = 0, 360, 5
	   xc = rxx + rrad * cos(ian * dtr)
	   yc = ryy + rrad * sin(ian * dtr)
	   val = integ(xc, yc, 4)
	   write(8, *) ian, val
	   if (val .gt. vmax) vmax = val
	   if (val .lt. vmin) vmin = val
	enddo
	type *, 'Min, max, %diff: ', vmin, vmax,
     *		200.0*(vmax-vmin)/(vmax+vmin)
c
	goto  1111
9999	end
c
c
c
	real*4 function integ( x, y, width)
	common ibuf
	integer*2 ibuf(2000,2000)
	real*4 x, y
	integer*4 width
c
	ix = (x + 0.5)
	iy = (y + 0.5)
	isum = 0
	icnt = 0
	do i1 = ix-width, ix+width
	   do i2 = iy-width, iy+width
	      isum = isum + ibuf(i1, i2)
	      icnt = icnt + 1
	   enddo
	enddo
	integ = float(isum) / icnt
	return
	end
c
c
c
	integer*4 function iexpand(i2)
	integer*2 i2
c
	if (i2 .eq. -1) then
	   type *, '***Overflow pixels (>65534)***'
	   stop '***Retake image with shorter exposure time'
	endif
	if (i2 .lt. 0) then
	   iexpand = 65536 + i2
	else
	   iexpand = i2
	endif
c
	return
	end