; Copyright (c) 1998-2006 A.P. Hitchcock  All rights reserved
;+
;NAME:
;	LITHO_READ
;
;LAST CHANGED: ----------------------------------- 22-feb-06
;
;PURPOSE:
;	This function reads a template (*.xyt)  image
; ( a set of (x,y,t) triplets ) for possible modification
;
;CATEGORY:
;	STAND ALONE: utility
;   called from aXis2000 as read~other~axis_xyt
;
;CALLING SEQUENCE:
;	Result = LITHO_READ(tmp [, NOF=nof, file=file])
;
;CALLED FROM AXIS:
;	Read->Images->AXIS
;
;INPUTS:
;	 ascii file with multi-column data
;
;KEYWORDS:
;	FILE     - filename to write
;   GROUP    - group of parent (axis_ID from axis)
;	VERBOSE  - print feedback on progress for trouble shooting
;
;OUTPUTS:
;	RESULT = TMP - axis format structure
;
;COMMON BLOCKS:
;	@AXIS_COM	standard set of common blocks
;
;;MODIFICATION HISTORY:
; (03-Nov-05 aph) developed from litho_save
; (22-feb-06 aph) fix array dimension erorr that prevented read of some xyt files (required kluge)
;-

function litho_read, file = file, group = group, verbose=verbose
@axis_com
on_error,2

axis_on=0
if NOT keyword_set(group) then group = 0 else begin
	t = size(axis_ID, /type)
	if t NE 0 then begin
		group = axis_ID
		axis_on=1
	endif else begin
		group = lox_ID
		axis_on=0
	endelse
endelse

; check if requested file exists; if not, get filename from user
IF keyword_set(file) then begin
	test = findfile(file)
  	if test(0) NE file then goto, get_filename
ENDIF ELSE BEGIN			; filename was not supplied or invalid
	get_filename:
	fltr='*.xyt'
	if keyword_set(nof) then fltr='*'
	file = pickfile2(/read, filter=fltr,  /LPATH, DEFPATH=defpath)
	if file EQ '' then return, file
ENDELSE

t = ax_name(file)
file = t(0)+t(1)+'.xyt'
fileshort = t(1)+'.xyt'

; ----- read the  (x,y,t) format file
openr, iunit, file, /get_lun

; -----------  read header ----------
line = ' '
while not  eof(iunit) do begin
	readf,iunit, line
	;print, line				; debug to ensure full file read-in
	test = strpos(line,'X-range')
	if test GT 0 then begin
		t_begin = strpos(line, '"', test-1)
		t_end = strpos(line, '"',t_begin+1)
		Xrange = float(strmid(line, t_begin+1, t_end-t_begin-1))
		test2 = strpos(line,'Y-range')
		t_begin = strpos(line, '"', test2-1)
		t_end = strpos(line, '"',t_begin+1)
		Yrange = float(strmid(line, t_begin+1, t_end-t_begin-1))
		test2 = strpos(line,'step')
		t_begin = strpos(line, '"', test2-1)
		t_end = strpos(line, '"',t_begin+1)
		step = float(strmid(line, t_begin+1, t_end-t_begin-1))
	endif
endwhile
close, iunit
free_lun, iunit

if keyword_set(verbose) then print,'X-range, Y-range, step ', Xrange, Yrange, step
WIDGET_CONTROL, /Hourglass
openr, iunit, file, /get_lun
; ----- skip to end of header ---------
while not  eof(iunit) do begin
	readf,iunit, line
	test = strpos(line,'/xyt-header')
	if test GT 0 then goto, end_header
endwhile
end_header:
j=long(0) & k=long(0)		; need long to avoid 32K limit
npts = 1
text = ''
x=fltarr(npts)
y=fltarr(npts)
t=fltarr(npts)
ix=lonarr(npts)				; need long to avoid 32K limit
iy=lonarr(npts)				; need long to avoid 32K limit
 while NOT eof(iunit) do begin
	readf, iunit,  text
	reads, text, xt, yt, tt, iyt, ixt   ; transpose needed since (iy, ix) listed in *.xyt
	x(j) = xt
	y(j) = yt
	t(j) = tt
	ix(j) = ixt
	iy(j) = iyt
	k = k+1
	if k EQ 10 then	begin
		if keyword_set(verbose) then print, j, x(j), y(j), t(j), ix(j), iy(j)
		k = 0
	endif
	j = j +1
	x = [x, 0.]
	y = [y, 0.]
	t=  [t, 0.]
	ix = [ix,0]
	iy = [iy,0]
;endfor
 endwhile
npts = j
close, iunit
free_lun, iunit

; ------ reconstruct original image based on file and user input
; xmin = min(x, max=xmax)
;
; ymin = min(y, max=ymax)
; npoints_x = max(ix)
; npoints_y = max(iy)
; npoints = max(npoints_x, npoints_y)

;if axis_on then xrange = get_num(prompt = "Xrange", val=range, group = axis_ID) $
;	 else xrange = get_num(prompt = "Xrange", val=xmin)
;if axis_on then yrangen = get_num(prompt = "Yrange", val=yrange, group = axis_ID) $
;	 else  yrange = get_num(prompt = "Yrange", val=yv)
;step = x(1) - x(0)
;if axis_on then step = get_num(prompt = "step size (nm)", val=step, group = axis_ID) $
;	 else  step = get_num(prompt = "step size (nm)", val=step)
;if axis_on then nx = get_num(prompt = "# of X-axis points (>)", val=nx, group = axis_ID) $
;	 else  ny = get_num(prompt = "# of Y-axis points (>)", val=ny)

; ---------------- set up arrays and build the image -------------------------
nx = fix(xrange/step)+3		; rounding error kluge !!
ny = fix(yrange/step)+3		; rounding error kluge !!
;print, 'npts, nx, ny', npts, nx, ny
;print, 'max (nx,ny)', max(nx), max(ny)
if  nx LT max(nx) then nx = max(nx)		; insurance !!
if  ny LT max(ny) then ny = max(ny)
tx = step*findgen(nx)
ty = step*findgen(ny)
d = fltarr(nx,ny)
for i = 0L, npts-1 do begin			; need long to avoid 32K limit
	d(ix(i),iy(i))= t(i)
endfor
energy = 0.0

tmp = {t:'2d', x:tx, y:ty, E: energy, d: d, xl:'  ', yl:'  ', dl: fileshort}
axis_log, 'read ' + string(npts) + ' points from litho file: ' + fileshort
axis_log, 'Full size of image is ' + string(tmp.x(n_elements(tmp.x)-1)) + $
            ' by '+ string(tmp.y(n_elements(tmp.y)-1))
return, tmp
END