       subroutine csrel(work,nwork,inpdev,lumdev,ierr)
	common /lumin/  sun(3), hdist, bmean, bsigma, bmult, boffst,
     1               zmult, zmean, zsigma
	common /bdist/  bmin, bmax, cntras,sunaz,sunel,linear,df(64)
c
       character*56    title
	character       pgm*8 
       dimension work(1)
c
       ierr=0
	bmult  = 0
	boffst = 0
       read(inpdev) title,pgm,nc,nr,nx,xo,dx,yo,dy
	if ((5* nc) .gt. nwork ) then
c  5*nc must be <= nwork for illuminated color slicing
         nw = nwork/5
	  print *, ' limited to ',nw,' columns'
         ierr=1
	  return
	endif
c
c  internal coordinates are x east, y north, z up.
	d2r    = 3.1415927 / 180.
	sun(1) = cos(d2r*sunel) * sin(d2r*sunaz)
	sun(2) = cos(d2r*sunel) * cos(d2r*sunaz)
	sun(3) = sin(d2r*sunel)
c
c  illumination
c  adjust z scaling from a sample of the input grid
	hdist = abs( 2. * dx )
	ipop  = 1 + (2*nc)
	npmax = 500
	call getpop( npop, npmax, nc, work(ipop), work, inpdev )
	ihist = ipop + 4 * npop
	sdeg = 7.
	call adjzm( sdeg, sunel, npop, work(ipop), work(ihist) )
	call adjdis( npop, work(ipop), work(ihist) )

c  generate illumination grid
	open( unit=lumdev, form='unformatted' )	
	it    = 1 + 3*nc
	call illum( nc, nr, work, work(it), inpdev, lumdev )

999	return
	end
c*******************************************************************************
	subroutine adjdis( npop, pop, brite )
	common /bdist/  bmin, bmax, cntras,sunaz,sunel,linear,df(64)
	dimension pop(4,npop), brite(npop)
c       character resp
c
c  initialize distribution function
1	call setdis
	do i = 1, npop
	  brite(i) = disfun( brite(i) )
	enddo
c
c  display illumination histogram
c      print *, 'illumination histogram'
c      call hist( npop, brite )
c
c      print *, ' present cntras =', cntras
c      write( *,'(a\)') ' do you want to test another ?[n]:'
c      read(*,'(a)') resp
c        if(resp.ne.'y'.and.resp.ne.'Y') go to 999
c      write( *,'(a\)') ' enter new cntras:'
c      read (*,*) cntras
c      call statis( npop, pop, brite )
c      go to 1
c
999	return
	end
 
c*******************************************************************************
	subroutine adjzm( sdeg, sunel, npop, t, g )
c  adjust scaling of z values to bring one standard deviation of the 
c  incident illumination angle distribution within sdeg.
c  
	common /lumin/ sun(3), hdist, bmean, bsigma, bmult, boffst, 
     1              zmult, zmean, zsigma
	dimension t(4,npop), g(npop)
	data d2r/ .01745329 /
c
	sd    = sin( d2r * sdeg ) * cos( d2r * .9 * sunel )
c
c  initial multiplier
	call ammi( 4*npop, t, tmin, tmax, mn, mx, 1 )
	tmax  = amax1( abs(tmax), abs(tmin) )
	zmin  = 0.
	zmult = 1.
	if ( tmax .gt. 10.*hdist .or. 
     1    tmax .lt. .01*hdist )     zmult = hdist / tmax
	zmax  = zmult
	call statis ( npop, t, g )
c
	if ( bsigma .lt. sd ) then
c  increase zmult until gradient distribution exceeds optimal
	do iter = 1, 10
	  zmult = 10. * zmult
	  call statis( npop, t, g )
	  if ( bsigma .gt. sd ) then
c  ok, we have an upper limit
	    zmax = zmult
	    go to 20
	    else
	    zmin = zmult
	  endif
	enddo
	endif
c
c  use bisection to refine zmult
20	do iter = 1, 10
	  zmult = .5 * ( zmin + zmax )
	  call statis( npop, t, g )
	  if ( bsigma .gt. sd ) then
	    zmax = zmult
	    else
	    zmin = zmult
	  endif
	enddo
c
	return
	end
 
c******************************************************************************
	subroutine ammi(n,a,amn,amx,mn,mx,isave)
	dimension a(n)
	if(isave.gt.1) go to 1
	mn=1
	mx=1
	amn=a(1)
	amx=a(1)
1	do 3 i=1,n
	if(a(i).ge.amn) go to 2
	amn=a(i)
	mn=i
2	if(a(i).le.amx) go to 3
	amx=a(i)
	mx=i
3	continue
	return
	end
 
 
c*******************************************************************************
	function disfun( x )
c  distribution function 
c  stretch brightness function to increase contrast
	common /lumin/  sun(3), hdist, bmean, bsigma, bmult, boffst,
     1               zmult, zmean, zsigma
	common /bdist/  bmin, bmax, cntras,sunaz,sunel,linear,df(64)
c
	disfun = bmult * ( x - boffst )
	if ( disfun .lt. bmin) then
	  disfun = bmin
	  else if ( disfun .gt. bmax ) then
	    disfun = bmax
	endif
c
	if ( linear .eq. 0 ) then
	  index = int ( 32. * ( disfun + 1.0 ) )
	  disfun = df( index )	
	endif
	return
	end	
 
c*******************************************************************************
	subroutine getpop( npop, ntmax, nc, t, z, isw )
c  sample input grid and store 4 element clusters
	dimension t(4,ntmax), z(nc)
	character id*56, pgm*8
	dval = 1.e38
	zchk = 1.e29
	do j = 1, ntmax
	  do i = 1, 4
	    t(i,j) = dval
	  enddo
	enddo
c
	rewind isw
	read(isw) id, pgm, nc1, nr, nz, xo, dx, yo, dy
c
	nskip = 1 + int( sqrt( float( nc*nr ) / float( ntmax ) ) )
10	nskip = nskip + 1
	nx = 1 + ( (nc-1) - 2 ) / nskip
	ny = 1 + ( (nr-1) - 2 ) / nskip
	if ( nx*ny .le. 1 ) stop ' nskip problem'
	if ( nx*ny .gt. ntmax ) go to 10
c
	if ( nskip .lt. 3 ) nskip = 3
	igap = nskip - 3
c
c  assemble population array  where:     t4
c                                      t1  t2
c                                        t3
c
	nt = 0
	do irow = 2, nr-1, nskip
c
	  call io( nc, z, -1, isw, isw, iend )
	  if ( iend .ne. 0 ) go to 50
	  nt2 = nt
	  do icol = 2, nc-1, nskip
	    nt2 = nt2 + 1
	    t(3,nt2) = z(icol)
	  enddo
c
	  call io( nc, z, -1, isw, isw, iend )
	  if ( iend .ne. 0 ) stop 111
	  nt2 = nt
	  do icol = 2, nc-1, nskip
	    nt2 = nt2 + 1
	    t(1,nt2) = z(icol-1)
	    t(2,nt2) = z(icol+1)
	  enddo
c
	  call io( nc, z, -1, isw, isw, iend )
	  if ( iend .ne. 0 ) stop 222
	  nt2 = nt
	  do icol = 2, nc-1, nskip
	    nt2 = nt2 + 1
	    t(4,nt2) = z(icol)
	  enddo
c
c  row of samples completed
	  nt = nt2
	  do i = 1, igap
	    read(isw, end=50)
	  enddo
	enddo
c
c  check each t group for dvals
50	jt = 0
	do it = 1, nt
	  iflag = 0
	  do i = 1, 4
	    if ( t(i,it) .gt. zchk ) iflag = 1
	  enddo
	  if ( iflag .eq. 0 ) then
c copy t(it) to output t(jt)
	    jt = jt + 1
	    do i = 1, 4
	      t(i,jt) = t(i,it)
	    enddo
	  endif
	enddo
c
	npop = jt
	if ( npop .le. 0 ) then
	  print *, ' nskip, raw npop, npopulation =', nskip, nt, npop
	  print *, ' statistical population empty'
	  print *, ' check your input grid for too many dvals.'
	  stop
	endif
	return
	end
 
c***********************************************************************
	subroutine io(n,z,iop,idev,jdev,iend)
c  read iop<0, write iop=0, r&w iop=1
	dimension z(n)
	iend=0
	if(iop)1,2,1
1	read(idev,end=10) y,z
	if(iop)9,9,2
2	write(jdev) y,z
9	return
10	iend=1
	return
	end
c*******************************************************************************
	subroutine hist( n, z )
	dimension zmm(2), ihist(20), z(n)
	character bar*66
	data nbmax, lenbar, dval  / 20, 66, 1.e38 /
	do i = 1, nbmax
	  ihist(i)=0
	enddo
	do i = 1, lenbar
	  bar(i:i)='*'
	enddo
	zmm(1) = dval
	zmm(2) = -zmm(1)
	do i = 1, n
	  if ( z(i) .lt. zmm(1) ) zmm(1) = z(i)
	  if ( z(i) .gt. zmm(2) ) zmm(2) = z(i)
	enddo
c
	call setdx( zmm, dz, nbmax )
	zmax = dz * aint( zmm(2) / dz + .49999 ) + dz / 2.
	zmin = dz * aint( zmm(1) / dz - .49999 ) - dz / 2.
	print *, 'zmin, zmax, dz=', zmin, zmax, dz
	nbar = ( zmax - zmin ) / dz  + 1
	if  ( nbar .gt. nbmax ) nbar = nbmax
c
	do i = 1, n
	  if ( z(i) .lt. dval ) then
	    k = 1 + int( ( zmax - z(i) ) / dz  + 0.5 )
	    if ( k .le. nbar .and. k .ge. 1 ) ihist(k) = ihist(k) + 1
	  endif
	enddo
c
	max=0
	do i = 1, nbar
	  if ( ihist(i) .gt. max ) max = ihist(i)
	enddo
c
	zlabel = zmax - dz / 2.
	do i = 1, nbar
	  nx = ( lenbar * ihist (i) ) / max
	  if ( ihist(i) .gt. 0 .and. nx .eq. 0 ) nx = 1
	  if ( nx .gt. 0 ) then
	    print 280, zlabel, ( bar(k:k), k = 1, nx )
280	    format( 1x, 1pg10.3, 2x, 66a1 )
	    else
	    print 280, zlabel
	  endif
	  zlabel = zlabel - dz
	enddo
c
	return
	end
 
c*******************************************************************************
	subroutine illum( nc, nr, z, t, isw, jsw )
c  dot the surface normal vector with the illumination vector.
c
	common /lumin/ sun(3), hdist, bmean, bsigma, bmult, boffst,
     1              zmult, zmean, zsigma
c
	dimension nn(3), xx(4), z(nc,3), t(nc)
	character id*56, pgm*8
	data fltbig/ 1.0e29 /, dval/ 1.0e38 /
	zfactr = zmult / hdist
c
	rewind isw
	read(isw) id, pgm, nn, xx
	write(id,5) sun, hdist
5	format( 'illum dircos=', 3f8.4, ' hd=', g15.5 )
	write(jsw) id, pgm, nn, xx
c
	do j = 1, 3
	  call io( nc, z(1,j), -1, isw, jsw, iend )
	  if ( iend .ne. 0 ) stop ' illum: eof in init buffer'
	enddo
c
	ifir = 1
	jrow = 0
	do irow = 2, nr-1
	  mid = ifir + 1
	  if ( mid .gt. 3 ) mid = 1
	  last = mid + 1
	  if ( last .gt. 3 ) last = 1
c
	  do 10 icol = 2, nc-1
	    t(icol) = dval
c
	         z1 = z(icol-1, mid)
	    if ( z1 .gt. fltbig ) go to 10
	         z2 = z(icol+1, mid)
	    if ( z2 .gt. fltbig ) go to 10
	         z3 = z(icol, ifir)
	    if ( z3 .gt. fltbig ) go to 10
	         z4 = z(icol, last)
	    if ( z4 .gt. fltbig ) go to 10
c
	    xdr = zfactr * ( z2 - z1 )
	    ydr = zfactr * ( z4 - z3 )
c	    xdr = zfactr * ( z(icol+1,mid) - z(icol-1,mid) )
c	    ydr = zfactr * ( z(icol, last) - z(icol,ifir)  )
	    fz  = 1.0 / sqrt( 1.0 + xdr*xdr + ydr*ydr )
	    fx  = -xdr * fz
	    fy  = -ydr * fz
	    traw = sun(1)*fx + sun(2)*fy + sun(3)*fz
	    t(icol) = disfun( traw )
10	  continue
	  t(1)  = t(2)
	  t(nc) = t(nc-1)
c
c  output illumination row
	  call io( nc, t, 0, isw, jsw, iend )
	  jrow = jrow + 1
         call perout(jrow,nr,1,1,'SHADED   ')
	  if ( irow .eq. 2 .or. irow .eq. nr-1 ) then
	    call io( nc, t, 0, isw, jsw, iend )
	    jrow = jrow + 1
	  endif
c
c  input next data row
	  if ( irow .lt. nr-1 ) then
	    call io( nc, z(1,ifir), -1, isw, jsw, iend )
	    if ( iend .ne. 0 ) print *, 'illum: eof at input row', irow
	    ifir = ifir + 1
	    if ( ifir .gt. 3 ) ifir = 1
	  endif
	enddo
c
	if ( jrow .ne. nr ) print *, 
     1    ' illum: input vs. output rows=', nr, jrow
       call perout(jrow,nr,1,1,'SHADED   ')
	return
	end
 
c*******************************************************************************
	subroutine setdis
c  set distribution function 
	common /lumin/  sun(3), hdist, bmean, bsigma, bmult, boffst,
     1               zmult, zmean, zsigma
	common /bdist/  bmin, bmax, cntras,sunaz,sunel,linear,df(64)
	linear = 1
c
	boffst = bmean
	bmult  = ( bmax - bmin ) / ( 4.0 * bsigma )
c
	if ( cntras .eq. 0 ) return
	linear = 0
	pi     = 3.1415927
c  modify a fraction of the distribution
	xd     = ( bmax - bmin ) / 2.0
	cd     = sqrt ( 2.0 * xd * xd )
	scaamp = ( cntras * cd ) / ( 2.0 * pi )
	cinc   = ( 2. * pi ) / ( 64. * xd / 2.0 ) 
	yinc   = 2. / 64.
c
	y   = -1.
	rad = -pi
	nrad  = int( 64. * xd / 2. )
	is    = ( 64 - nrad ) / 2
	ie    = is + nrad
c
	do i = 1, 64
	  df(i) = y
	  if ( i .ge. is .and. i .le. ie ) then
	    c     = scaamp * sin( rad ) 
	    df(i) = df(i) + c
	    rad   = rad + cinc
	  endif
	  y = y + yinc
	enddo
	return
	end	
 
c*******************************************************************************
	subroutine setdx( range, dx, maxint )
	dimension range(2)
	ixpon(r) = int( alog10( abs( r ) ) + 100. ) - 100
	maxi     = maxint
	if ( maxi .le. 0 ) maxi = 20
	t   = abs( ( range(2) - range(1) ) / float( maxi ) )
	if ( t .lt. 1.e-20 ) t = 1.e-20
	p10 = 10.0 ** ixpon( t )
	t1  = t / p10
	if ( t1 .le. 1.0 ) dx =       p10
	if ( t1 .gt. 1.0 ) dx =  2. * p10
	if ( t1 .gt. 2.0 ) dx =  5. * p10
	if ( t1 .gt. 5.0 ) dx = 10. * p10
	if ( t  .le. 1.e-20 ) dx = 1.
	if ( range(1) .gt. range(2) ) dx = -dx
	return
	end
 
c*******************************************************************************
	subroutine statis( npop, t, g )
	common /lumin/ sun(3), hdist, bmean, bsigma, bmult, boffst,
     1              zmult, zmean, zsigma
	dimension t(4,npop), g(npop)
	double precision gs, gsq
	if ( zmult .eq. 0.0 ) zmult = 1.0
c
	zfactr = zmult / hdist
	gs  = 0.d0
	gsq = 0.d0
	do j = 1, npop
	  xg = zfactr * ( t(2,j) - t(1,j) )
	  yg = zfactr * ( t(4,j) - t(3,j) )
	  cosz = 1.0 / sqrt( 1 + xg*xg + yg*yg )
	  cosx = -xg * cosz
	  cosy = -yg * cosz
	  g(j) = sun(1)*cosx + sun(2)*cosy + sun(3)*cosz
	  gs   = gs  + dble( g(j) )
	  gsq  = gsq + dble( g(j) ) * dble( g(j) )
	enddo
c
	gpop   = npop
	sn     = sngl( dble(gpop) * gsq - gs * gs ) 
	bsigma = 0.0
	if ( sn .gt. 0.0 ) bsigma = sqrt( sn / (gpop * (gpop - 1.0)) )
	bmean  = sngl(gs) / gpop
	return
	end
