;file processmovies.ana, ras
; /3/93 second generation movie processer adapted from alignmovies
 ;user interface to generate parameters
 ;
 ;9/23/90	rotation upgraded and reference angle allowed, if r_ref
 ;		is defined in command file, it is used rather than computing
 ;		a mid range, also lt 0 forces compute -- in a series the
 ;		first movie will define r_ref and the other will use the same
 ;5/23/90	some changes to handle 1024 images, these involve assuming
 ;that images are 512 until the first is read in, then we change the
;subarea (for CC) and the rotation center (if a rotation is done)
 ;
 ;an intermediate step to a better system, this file contains routines to
 ;rigidly align a series of images (stored as disk files in FZ format) and
 ;to destretch them
 ;the routines are run via command files which define the file names,
 ;length of sequence, and operation performed
 ;
 ;both the alignment and the destretch use a double pass method, the first
 ;pass determines the image by image offsets and the second does some trend
 ;processing on these before applying them and producing output images
 ;file name conventions are assumed to comply with the following rules:
 ;	last 5 characters of file name are n#### were the #'s represent the
 ;	number of the image in the sequence; i.e., the first file name
 ;	might be LB113w5576N0001. The first part of the file name is
 ;	called name_set within the software. The files must progress from
 ;	0001 up to the number of files in steps of 1.
 ;	The file extensions are .COR for the initial images. The software
 ;	then creates .RIG files and .STRT files for the rigidly aligned and
 ;	destretched results.
@/umbra/people/shine/auto/series_align_new
@/umbra/people/shine/auto/rigplot
@/umbra/people/zoe/look92
 ;==============================================================================
block getstarted
 close,1,2
 #diskout1='/data2/soup/'
 #diskout2='/data0/soup/'
 #diskout3='/data1/soup/'
 #diskout4='/pore2/people/shine/lp92/'
 #diskout5='/data5/soup/'
 #diskout6='/pore2data/soup/'
 #diskout7='/data/soup/'
 #ndisksout=7		;number to use
 ;get id for this data set
 tapeid=''
 read,'enter catalog tape id (e.g., le132f18)',tapeid
 catload,tapeid
 type,'nimage =',nimage,', nexp =',nexp,'
 if nexp gt 0 then { nrepeats=nimage/nexp
	  type,'nrepeats (complete ones) =',nrepeats }
 ;check for duplicates that should be combined
 waveindx=lonarr(nexp)
 nc=1
 ;loop through and assign a number for each image, images with the same state
 ;get the same number
 for k=0,nexp-1 do {
 ;gen the name, assume unique and then check
 ss='namew'+istring(nc,3,2)
 wname,seq_wave(k),seq_off(k),seq_pol(*,k),eval(ss),tapeid
 ;now check against all previous names
 waveindx(k)=nc	nc=nc+1
 if k gt 0 then { for i=0,k-1 do { s2='namew'+istring(waveindx(i),3,2)
   if eval(s2) eq eval(ss) then { nc=nc-1  waveindx(k)=waveindx(i)  break }}}
 }
 nwave=nc-1
 ty,'number of unique wavelengths is',nwave
 ;get subcycles for each wavelength
 subcycle=zero(intarr(nwave))
 for k=0,nexp-1 do {
 subcycle(waveindx(k)-1) += 1
 }
 endblock
 ;==============================================================================
block autotwist 
 ty,'please select an option'
 ty,'0: compute alignments and produce output files'
 ty,'1: only compute alignments, defer output'
 ty,'2: only produce output files (assumes alignments on disk)'
 ty,'3: rotate but don''t align, produces output files'
 read,#twist_option
 if #twist_option eq 0 then ty,'if things get dicey, output might be deferred'
 run getstarted		;load the catalog info and such
 
 wavesave=zero(intarr(nwave))+1
 mess=''
 read,'are we processing all the wavelengths ? [Y/N]',mess
 if upcase(mess) ne 'Y' then {
 type,'for each wavelength set, type Y to untwist, N to skip, Q to skip all remaining'
 zero,wavesave		;a flag array
 for iw=1,nwave do {
 type,eval('namew'+ist(iw,3,2))
 read,mess
 if upcase(mess) eq 'Y' then wavesave(iw-1)=1
 if upcase(mess) eq 'Q' then break
 }
 }
 kfirstrep=1
 kreps=nrepeats
 while 1 eq 1 do {
 mess=''
 read,'Are we processing all the repeats ? [Y/N]',mess
 if upcase(mess) eq 'Y' then break
 read,'enter first and last rep (of the complete cycles) to use',kfirstrep,kreps
 type,'you have first and last rep =', kfirstrep,kreps
 read,'is this OK [Y/N]?',mess
 if upcase(mess) eq 'Y' then break
 type,'let''s try again then'
 }
 repeats_used=kreps-kfirstrep+1
	 ;set up some parameters
 #rotate=1		;we are rotating by default (chance to change later)
 r_ref=-1.0		;use mean rotation by default
	 ;determine drift parameters, first get duration of run
 dt=3600.*fix(#t(39:40))+60.*fix(#t(42:43))+fix(#t(45:46)
 dt=dt-3600.*fix(#t(27:28))-60.*fix(#t(30:31))-fix(#t(33:34)
 if dt lt 0 then dt=dt+24.*3600.
 dt_rep=dt/nrepeats
	 ;if only a partial run, adjust dt
 dt=dt_rep*(repeats_used)
   ;our scheme - if the duration < 30m, don't use superdrift, just rigid align
   ;if > 30 but < 90, smooth with about a 10 minute window, for longer
   ;runs use 15 m, these may need adjusting as experience improves
   ;also use rigid (no drift) if repeats_used is 10 or less
 if dt le 1800 or repeats_used le 10 then { driftbase=0 } else {
 if dt le 5400 then { driftbase=600./dt_rep } else { driftbase=900./dt_rep }
 }
 ty,'computed driftbase =', driftbase
	 ;choose an output disk
 disk=zero(lonarr(nwave))	;provide for different disks later
 
 if #twist_option ne 1 then {
 iq=-1
 while (iq le 0 or iq gt #ndisksout) do {
 ty,'enter a number to choose output disk'
 for i=1,#ndisksout do { ty,i,'  ',eval('#diskout'+ist(i,2,2)) }
 read,iq
 }
 disk=disk+iq-1
 
 read,'enter 1 to turn on delete flag', #dflag
 if #dflag ne 1 then #dflag=0
 }
 
 name_in_ext='.cor'
 name_out_ext='.rig'
 mess=''
 ty,'anything special ?'
 ty,'(turn off rotation, override drift factor, non-standard name extensions,'
 ty,'special subarea)'
 read,'[Y/N]',mess
 if upcase(mess) eq 'Y' then {
 read,'enter 1 for La Palma rotation or 0 for none', #rotate
 #rotate=fix(#rotate)
 mess=''
 if #rotate ne 0 then {
 read,'enter rotation reference angle or -1 to use mean orientation',r_ref
 }
 read,'want to change driftbase ? [Y/N]',mess
	 if upcase(mess) eq 'Y' then {
	 ty,'enter new driftbase (0 for completely rigid)', driftbase }
 read,'want to change input name extension ? [Y/N]',mess
	 if upcase(mess) eq 'Y' then {
	 read,'enter it with the "." (e.g., ".cor")', name_in_ext }
 read,'need to specify the CC subarea ? [Y/N]',mess
	 if upcase(mess) eq 'Y' then {
	 ty,'please be very careful !'
	 read,'enter first and last x index',#i1,#i2
	 read,'enter first and last y index',#j1,#j2
	 #usersubarea=1
	 }
 if #twist_option ne 1 then {
 read,'want to change output name extension ? [Y/N]',mess
	 if upcase(mess) eq 'Y' then {
	 read,'enter it with the "." (e.g., ".rig")', name_out_ext }
 }
 
 }
 
 ;log our beginning
 mess=istring(fix(!systime),10,2)	;a unique part of the file name
 close,1	openw,1,lowcase(tapeid)+'.'+mess+'.history'
 printf,1,tapeid,' de-rotation and de-jitter processing begun ',!date,' ',!time
 ;also print info block on plot file
 pdev,1
 run infoblock
 
 count=0
 ;main untwist loop
 
 for iw=0,nwave-1 do {
 if wavesave(iw) ne 0 {
 
   if driftbase gt 0 then {
   #drift=fix(driftbase*subcycle(iw))
   #drift=#drift-(#drift%2 eq 0)			;make odd
   #drift=#drift>3
   } else #drift=0
 
   count = count + 1
   ;get the file counts
   kstart=1+ (kfirstrep-1)*subcycle(iw)
   kend=kreps*subcycle(iw)
   name_set=lowcase(eval('namew'+ist(iw+1,3,2)))
   name_in=name_set+'n####'+name_in_ext
   printf,1,!time
   printf,1,'iw =', iw,' kstart and kend =', kstart,kend
   printf,1,'name_in =', name_in
   ty,'iw =', iw,' kstart and kend =', kstart,kend
   ty,'name_in =', name_in

 ;need angles for all cases but 2 (which reads in old values)
 if #twist_option ne 2 then {
   getangles,kstart,kend,name_in,r_ref,r
 }
 
 ;options, getrigid, saves, and plots done only for 0 and 1
 if #twist_option eq 0 or #twist_option eq 1 then {
 
   getrigid,kstart,kend,name_in,r_ref,r,d,seeing,meanint,sample_a,sample_b
 
	 ;getrigid computes the angles and the pair offsets
	 ;saved in r and d
	 ;note that #twist_option can be changed to 1 by getrigid
 
   s='angles from autoprocess, first value is ref, remainder has it subtracted'
   fzwrite,[r_ref,r],name_set+'.angles',s
   s='offsets from getrigid in autoprocess'
   fzwrite,d,name_set+'.offsets',s
   s='seeing measure from getrigid in autoprocess'
   fzwrite,seeing,name_set+'.seeing',s
   ;plot some results
   if count eq 1 then  { rotplot,r_ref,r,name_set,kstart }
   ;also postage stamps of first and last images
   stamps,sample_a,sample_b,count,1,name_in_ext,kstart,kend
   ;before plotting, process the raw offsets
   dhipass,d,dxs,dys
   rigplot,count,name_set,dxs,dys,seeing,meanint,kstart
 } 
 ;for option 2, the offsets and angles are restored
 if #twist_option eq 2 then {
   fzread,r,name_set+'.angles'
   r_ref=r(0)		r=r(1:*)
   fzread,d,name_set+'.offsets'
   ;process the raw offsets
   dhipass,d,dxs,dys
 }

 ;for option 3, we make a zero delta array
 if #twist_option eq 3 then {
 nf=kend-kstart+1
 dxs=zero(fltarr(nf))		dys=dxs
 }
 
   sdisk=eval('#diskout'+ist(disk(iw)+1,2,2))
 ;this info printed even if we aren't writing files
   name_out=sdisk+name_set+'n####'+name_out_ext
   printf,1,'#twist_option =', #twist_option
   printf,1,'name_out =', name_out
   printf,1,'drift factor =',#drift
   ty,'#twist_option =', #twist_option
   ty,'name_out =', name_out
 
 ;note that untwist produces the samples for #twist_option = 1
 ;these are not printed for #twist_option = 2 or 3
 
 untwist,kstart,kend,name_in,name_out,r_ref,r,dxs,dys,sample_a,sample_b
 
 if #twist_option eq 0 or #twist_option eq 1 then {
   stamps,sample_a,sample_b,count,0,name_out_ext,kstart,kend
   ;need to eject?
   if count%3 eq 2 then erase
 }
 }}	;end of iw loop
 ;change name of plot file
 spawn,'mv junk.eps '+lowcase(tapeid)+'.twist.eps'
 printf,1,'untwist processing completed ',!date,' ',!time
 close,1
 endblock
 ;==============================================================================
block customtwist
 name_set=''	name_in_ext=''	name_out_ext=''
 #samples = 0	;causes no plots or little pictures
 read,'enter set name (without extension, e.g., le132f18w6302m350n####)',name_set
 read,'enter input extension (with the "." please)',name_in_ext
 read,'enter output extension (with the "." please)',name_out_ext
 read,'enter first and last file numbers', kstart,kend
 nk=kend-kstart+1
 if nk le 0 then {ty,'last must be ge than first' retall }
 read,'enter drift smoothing factor', #drift
 read,'enter 1 for La Palma rotation or 0 for none', #rotate
 #rotate=fix(#rotate)
 mess=''
 if #rotate ne 0 then {
 read,'enter rotation reference angle or -1 to use mean orientation',r_ref
 }
 
 read,'enter 1 to turn on delete flag', #dflag
 if #dflag ne 1 then #dflag=0
 ty,'please select an option'
 ty,'0: compute alignments and produce output files'
 ty,'1: only compute alignments, defer output'
 ty,'2: only produce output files (assumes alignments on disk)'
 ty,'3: rotate but don''t align, produces output files'
 read,#twist_option
 if #twist_option eq 0 then ty,'if things get dicey, output might be deferred'
 mess=''
 read,'need to specify the rotation center ? [Y/N]',mess
	 if upcase(mess) eq 'Y' then {
	 ty,'please be very careful !'
	 read,'enter rotation point x and y',#xcenter,#ycenter
	 #xcenter=#xcenter/2	#ycenter=#ycenter/2
	 #usercenter=1
	 }
 read,'need to specify the CC subarea ? [Y/N]',mess
	 if upcase(mess) eq 'Y' then {
	 ty,'please be very careful !'
	 read,'enter first and last x index',#i1,#i2
	 read,'enter first and last y index',#j1,#j2
	 #i1=#i1/2	#i2=#i2/2	#j1=#j1/2	#j2=#j2/2
	 #usersubarea=1
	 }
 
 ;log our beginning
 mess=istring(fix(!systime),10,2)	;a unique part of the file name
 close,1
 sq='custom'
 openw,1,lowcase(sq)+'.'+mess+'.history'
 printf,1,sq,' de-rotation and de-jitter processing begun ',!date,' ',!time
 
    name_in=name_set+name_in_ext
   printf,1,!time
   printf,1,' kstart and kend =', kstart,kend
   printf,1,'name_in =', name_in
   ty,' kstart and kend =', kstart,kend
   ty,'name_in =', name_in

 ;need angles for all cases but 2 (which reads in old values)
 ;unless no rotation of course
 if #rotate ne 0 then {
 if #twist_option ne 2 then {
   getangles,kstart,kend,name_in,r_ref,r
 } } else r=zero(fltarr(nk))
 
 ;options, getrigid, saves, and plots done only for 0 and 1
 if #twist_option eq 0 or #twist_option eq 1 then {
 
   getrigid,kstart,kend,name_in,r_ref,r,d,seeing,meanint,sample_a,sample_b
 
	 ;getrigid computes the angles and the pair offsets
	 ;saved in r and d
	 ;note that #twist_option can be changed to 1 by getrigid
 
   s='angles from autoprocess, first value is ref, remainder has it subtracted'
   fzwrite,[r_ref,r],lowcase(strr(name_set,'#',''))+'.angles',s
   s='offsets from getrigid in autoprocess'
   fzwrite,d,lowcase(strr(name_set,'#',''))+'.offsets',s
   s='seeing measure from getrigid in autoprocess'
   fzwrite,seeing,lowcase(strr(name_set,'#',''))+'.seeing',s
   ;before plotting, process the raw offsets
   dhipass,d,dxs,dys
 } 
 ;for option 2, the offsets and angles are restored
 if #twist_option eq 2 then {
   if #rotate ne 0 then {
    fzread,r,lowcase(strr(name_set,'#',''))+'.angles' } else r=zero(fltarr(nk))
   r_ref=r(0)		r=r(1:*)
   fzread,d,lowcase(strr(name_set,'#',''))+'.offsets'
   ;process the raw offsets
   dhipass,d,dxs,dys
 }

 ;for option 3, we make a zero delta array
 if #twist_option eq 3 then {
 nf=kend-kstart+1
 dxs=zero(fltarr(nf))		dys=dxs
 }
 
 ;this info printed even if we aren't writing files
   name_out=name_set+name_out_ext
   printf,1,'#twist_option =', #twist_option
   printf,1,'name_out =', name_out
   printf,1,'drift factor =',#drift
   ty,'#twist_option =', #twist_option
   ty,'name_out =', name_out
 
 ;note that untwist produces the samples for #twist_option = 1
 ;these are not printed for #twist_option = 2 or 3
 
 if #twist_option ne 1 then {
 untwist,kstart,kend,name_in,name_out,r_ref,r,dxs,dys,sample_a,sample_b
 }
 
 printf,1,'untwist processing completed ',!date,' ',!time
 close,1
 endblock
 ;=============================================================================
block autostretch
 ty,'please select an option'
 ty,'0: compute distortion grids and produce output files'
 ty,'1: only compute distortion grids, defer output'
 ty,'2: only produce output files (assumes distortion grids on disk)'
 read,#stretch_option
 run getstarted		;load the catalog info and such
 
 wavesave=zero(intarr(nwave))
 
 ;for destretching, always query for which wavelengths since we rarely
 ;do all of them
 
 type,'for each wavelength set, type Y to destretch, N to skip, Q to skip all remaining'
 for iw=1,nwave do {
 type,eval('namew'+ist(iw,3,2))
 mess=''
 read,mess
 if upcase(mess) eq 'Y' then wavesave(iw-1)=1
 if upcase(mess) eq 'Q' then break
 }
 
 kfirstrep=1
 kreps=nrepeats
 while 1 eq 1 do {
 mess=''
 read,'Are we processing all the repeats ? [Y/N]',mess
 if upcase(mess) eq 'Y' then break
 read,'enter first and last rep (of the complete cycles) to use',kfirstrep,kreps
 type,'you have first and last rep =', kfirstrep,kreps
 read,'is this OK [Y/N]?',mess
 if upcase(mess) eq 'Y' then break
 type,'let''s try again then'
 }
 repeats_used=kreps-kfirstrep+1
 
	 ;set up some parameters
	 ;determine drift parameters, first get duration of run
 
 dt=3600.*fix(#t(39:40))+60.*fix(#t(42:43))+fix(#t(45:46)
 dt=dt-3600.*fix(#t(27:28))-60.*fix(#t(30:31))-fix(#t(33:34)
 if dt lt 0 then dt=dt+24.*3600.
 dt_rep=dt/nrepeats
	 ;if only a partial run, adjust dt
 dt=dt_rep*(repeats_used)
   ;computed differently for stretching than for aligning
   ;since we may want to preserve flows even in short runs
   ;smoothing is done with a gaussian kernel, upgrade from boxcar
   ;our scheme - if the duration < 20m, use half the duration
   ;if > 20 but < 90, smooth with about a 10 minute window, for longer
   ;runs use 15 m, these may need adjusting as experience improves
   ;also use rigid (no drift) if repeats_used is 10 or less
 if dt le 1200 or repeats_used le 10 then { driftbase=repeats_used/2 } else {
 if dt le 5400 then { driftbase=600./dt_rep } else { driftbase=900./dt_rep }
 }
 ty,'computed driftbase =', driftbase
	 ;choose an output disk
 disk=zero(lonarr(nwave))	;provide for different disks later
 if #stretch_option eq 0 or #stretch_option eq 2 then {
 iq=-1
 while (iq le 0 or iq gt #ndisksout) do {
 ty,'enter a number to choose output disk'
 for i=1,#ndisksout do { ty,i,'  ',eval('#diskout'+ist(i,2,2)) }
 read,iq
 }
 disk=disk+iq-1
 read,'enter 1 to turn on delete flag', #dflag
 if #dflag ne 1 then #dflag=0
 }
 
 if #stretch_option eq 0 or #stretch_option eq 1 then {
 nest = 3
 grid_set=[8,16,32]
 grid_clips=[5,4,3]
 ty,'the default destretching is:  nest=',nest
 ty,'grid_set =',grid_set
 ty,'grid_clips =',grid_clips
 mess=''
 read,'enter 1 to change any of these',mess
 if mess eq '1' then {
  read,'enter number of nested destretches', nest
  grid_set=fltarr(nest)
  grid_clips=fltarr(nest)
  read,'enter grid sizes', grid_set
  read,'enter grid clips', grid_clips
 }
 } 

 read,'enter 1 to trim edges, 0 to keep the whole image', #trim_edge
 if #trim_edge ne 1 then #trim_edge=0
 
 name_in_ext='.rig'
 name_out_ext='.strt'
 mess=''
 ty,'anything special ?'
 ty,'(override drift factor, non-standard name extensions)'
 read,'[Y/N]',mess
 if upcase(mess) eq 'Y' then {
 mess=''
 read,'want to change driftbase ? [Y/N]',mess
	 if upcase(mess) eq 'Y' then {
	 read,'enter new driftbase (0 for completely rigid)', driftbase }
 read,'want to change input name extension ? [Y/N]',mess
	 if upcase(mess) eq 'Y' then {
	 read,'enter it with the "." (e.g., ".rig")', name_in_ext }
 read,'want to change output name extension ? [Y/N]',mess
	 if upcase(mess) eq 'Y' then {
	 read,'enter it with the "." (e.g., ".strt")', name_out_ext }
 }
 ;log our beginning
 mess=istring(fix(!systime),10,2)	;a unique part of the file name
 close,1	openw,1,lowcase(tapeid)+'.'+mess+'.history'
 printf,1,tapeid,' de-stretch processing begun ',!date,' ',!time
 !match_messages=0
 
 ;main loop
 for iw=0,nwave-1 do {
 if wavesave(iw) ne 0 {
   ;get the file counts
   kstart=1+ (kfirstrep-1)*subcycle(iw)
   kend=kreps*subcycle(iw)
   name_in=eval('namew'+ist(iw+1,3,2))+'n####'+name_in_ext
   printf,1,!time
   printf,1,'iw =', iw,' kstart and kend =', kstart,kend
   printf,1,'name_in =', name_in
   ty,'iw =', iw,' kstart and kend =', kstart,kend
   ty,'name_in =', name_in
  if #stretch_option eq 0 or #stretch_option eq 1 then {
   getstretch,kstart,kend,name_in,delta,grid_set,grid_clips,seeing
 
   ;save some results
   s='grid from getstretch in autoprocess'
   name_out=eval('namew'+ist(iw+1,3,2))
   fzwrite,delta,lowcase(name_out)+'.grid',s
  }

  if #stretch_option eq 2 then {
   ;read the grids back in
   name_out=eval('namew'+ist(iw+1,3,2))
   fzread,delta,lowcase(name_out)+'.grid',s
  }

   if driftbase gt 0 then {
   #drift=fix(driftbase*subcycle(iw))
   #drift=#drift-(#drift%2 eq 0)		;make odd
   #drift=#drift>3
   } else #drift=0

  if #stretch_option eq 0 or #stretch_option eq 2 then {
   sdisk=eval('#diskout'+ist(disk(iw)+1,2,2))
   name_out=sdisk+eval('namew'+ist(iw+1,3,2))+'n####'+name_out_ext
   printf,1,'name_out =', name_out
   printf,1,'drift factor =',#drift
   ty,'name_out =', name_out
   destretch,kstart,kend,name_in,name_out,delta
  }
 
 }}	;end of iw loop
 spawn,'mv junk.eps '+tapeid+'.stretch.eps'
 printf,1,'destretch processing completed ',!date,' ',!time
 close,1
 endblock
 ;==============================================================================
block customstretch
 close,1,2
 ty,'please select an option'
 ty,'0: compute distortion grids and produce output files'
 ty,'1: only compute distortion grids, defer output'
 ty,'2: only produce output files (assumes distortion grids on disk)'
 read,#stretch_option
 name_set=''	name_in_ext=''	name_out_ext=''
 read,'enter set name (without extension, e.g., le132f18w6302m350n####)',name_set
 read,'enter input extension (with the "." please)',name_in_ext
 read,'enter output extension (with the "." please)',name_out_ext
 read,'enter first and last file numbers', kstart,kend
 if #stretch_option eq 0 or #stretch_option eq 2 then {
 read,'enter smoothing factor', #drift
 }
 if #stretch_option eq 0 or #stretch_option eq 1 then {
 read,'enter number of nested destretches', nest
 grid_set=fltarr(nest)
 grid_clips=fltarr(nest)
 read,'enter grid sizes', grid_set
 read,'enter grid clips', grid_clips
 read,'enter maximum total shift', #max_stretch_shift
 }
 read,'enter 1 to trim edges, 0 to keep the whole image', #trim_edge
 if #trim_edge ne 1 then #trim_edge=0
 if #stretch_option eq 0 or #stretch_option eq 2 then {
	 ;choose an output disk
 #diskout1='/data2/soup/'
 #diskout2='/data0/soup/'
 #diskout3='/data1/soup/'
 #diskout4='/pore2/people/shine/lp92/'
 #diskout5='/data5/soup/'
 #diskout6='/pore2data/soup/'
 #diskout7='/data/soup/'
 #diskout8='/pore2data/'
 #diskout9='/pore2data2/'
 #ndisksout=9		;number to use
 iq=-1
 while (iq le 0 or iq gt #ndisksout) do {
 ty,'enter a number to choose output disk'
 for i=1,#ndisksout do { ty,i,'  ',eval('#diskout'+ist(i,2,2)) }
 read,iq
 }
 disk=iq
 read,'enter 1 to turn on delete flag', #dflag
 if #dflag ne 1 then #dflag=0
 }
 ;log our beginning
 mess=istring(fix(!systime),10,2)	;a unique part of the file name
 close,1
 openw,1,lowcase(strr(name_set,'#',''))+'.'+mess+'.history'
 printf,1,name_set,' de-stretch processing begun ',!date,' ',!time
 name_in=name_set+name_in_ext
  if #stretch_option ne 2 then {
  printf,1,'nest factor =',nest
  printf,1,'grid sizes =', grid_set
  printf,1,'grid clips =', grid_clips }
  if #stretch_option eq 0 or #stretch_option eq 1 then {
   ty,'name_in =', name_in, ' looping from',kstart,' to', kend
   printf,1,'name_in =', name_in
   getstretch,kstart,kend,name_in,delta,grid_set,grid_clips
 
   ;save some results
   s='grid from getstretch in autoprocess'
   ;special for running multiple sets, add time to file name
    mess=istring(fix(!systime),10,2)	;a unique part of the file name

   fzwrite,delta,lowcase(strr(name_set,'#',''))+'.grid'+mess,s
 }
  if #stretch_option eq 2 then {
   ;read the grids back in
   fzread,delta,lowcase(strr(name_set,'#',''))+'.grid',s
  }
  if #stretch_option eq 0 or #stretch_option eq 2 then {
   sdisk=eval('#diskout'+ist(disk,2,2))
   name_out=sdisk+name_set+name_out_ext
   printf,1,'name_out =', name_out
   printf,1,'drift factor =',#drift
   ty,'name_out =', name_out
 
   destretch,kstart,kend,name_in,name_out,delta
 } 
 printf,1,'destretch processing completed ',!date,' ',!time
 close,1
 endblock
 ;=============================================================================
subr wname,wave,offset,pol,name,tapeid
 ;does the wavelength, etc portion of name
 ;wave and offset are binary and pol is a string or byte array of length 3
 s='W'+istring(wave,4,2)
 case
 offset lt 0:	s=s+'M'+istring(abs(offset),4,2)
 offset eq 0:	s=s+'LC'
 offset gt 0:	s=s+'P'+istring(abs(offset),4,2)
 endcase
 sp=bmap(pol)	sp=sieve(sp, sp ne 32)	;remove blanks
 sp=smap(sp)
 if sp ne '0' then s=s+sp
 name=tapeid+s
 name=lowcase(name)
 endsubr
 ;==============================================================================
subr processmoviesabort
 printf,1,'processing aborted ',!date,' ',!time
 close,1
 retall
 endsubr
 ;==============================================================================
subr getangles(kstart,kend,name_in,r_ref,r)
 ;for VTT data and already derotated movies, disable rotation
 nf=kend-kstart+1
 if #rotate ne 0 then {
 rag=0
 lpgetangles,rag,name_in,kstart,kend	;get angles
 rag=rag*#r.d	;the rotate function just changes it back to radians of course
 if r_ref lt 0.0 then {		;9/23/90 - allow choice, figure r_ref only
				 ;if it has been set to lt 0
 r_ref=0.5*(max(rag)+min(rag)) }	;use the mid point for ref, another
 			 ;approach might be to use the time of best seeing
 r=rag-r_ref
 ;special for george's data, the image is reversed (I guess) and the angles
 ;are negative
 r=-r
 t,'reference angle =',r_ref
 } else { r=zero(fltarr(nf)) }
 endsubr
 ;==============================================================================
subr getrigid(kstart,kend,name_in,r_ref,r,deltac,seeing,meanint,s1,s2)
 ;2/3/93 subroutine version adapted from block version in alignmovies
 ;it does a NN rotate to save time and sub-pixel rigid align, saves alignment
 ;results for use in untwist
 
 if defined(#usersubarea) eq 0 then #usersubarea=0
 nf=kend-kstart+1
 deltac=zero(fltarr(2,nf))
 ;2/9/93 keep large shifts in a separate array
 shiftc=zero(fltarr(2,nf))
 seeing=fltarr(nf)
 meanint=fltarr(nf)
 lowcut=2	;removes the lower 2 freqs from the cross correlation
 hicut=0		;if hicut=0, then none applied
 m2=0		head=0		x=0		m1=0
 dxp=0.0		dyp=0.0
 ;loop over the images now
 
 ;-----------------------------------------------------------------------------
 for k=kstart,kend do {
 
 kmo=k-kstart
 fseek,xfull,fns(name_in,k),head		;read in the image
 x = compress(xfull,2)				;compress right away
 ;work only with x2 compressed all the way through since we use it for CC
 ;get sample images for kstart and kend
 
 if #samples eq 1 then {
 	if k eq kstart then s1=ft(compress(x,2))
 	if k eq kend then s2=ft(compress(x,2))
 }
 
 ;for the first value of k, set up subarea parameters
 if k eq kstart then { 
   nx=dimen(x,0)  ny=dimen(x,1)
   if #usercenter eq 0 then {	xrc=nx/2	yrc=ny/2 } else {
   				xrc = #xcenter	yrc = #ycenter  }
   if #usersubarea eq 0 then {	;default choices if not user specified
   nq=nx<ny  nsub=384
   if nq gt 640 then {	nsub=640 }
   ;other possibilities considered are 256, 196, 160, and 128 for small arrays
   if nq lt 384 then { xq=[64,92,128,160,196,256] 
	 nsub=max(sieve(xq, xq le nq))
	 if nsub le 0 then { t,'arrays too small for ALIGNMOVIES' retall } }
   t,'selected subarea size =',nsub
   ;use a centered subarea, normal is 384x384 or 640x640
   i1=(nx-nsub)/2	i2=i1+nsub-1	j1=(ny-nsub)/2	j2=j1+nsub-1
   } else {  i1=#i1  i2=#i2  j1=#j1  j2=#j2
   }	;end of #usersubarea conditional
   t,'selected subarea ranges are',i1,i2,j1,j2,' center',xrc,yrc
   printf,1,'selected subarea ranges are',i1,i2,j1,j2,' center',xrc,yrc	}
  ;note that center and subarea indices are for x2 compressed
 ;2/17/93 - seeing measure added, not very efficient, consider using
 ;something from transform space since we do an fft anyhow
 
 if #samples eq 1 then {
 xq=x(i1:i2,j1:j2)
 meanint(kmo)=mean(xq)			;mean intensity
 seeing(kmo)=mean(abs(differ(xq)))/meanint(kmo)	;"seeing" measure
 }
 
 ;rotate ?
 if #rotate ne 0 then m2=rotate(x,r(kmo),xrc,yrc) else switch,m2,x
 ;if this is not the first, do the alignment
 if k eq kstart then {
 ;for the first one we need to set up alignment code
 start_align,m2(i1:i2,j1:j2),lowcut,hicut
 ix=0	iy=0
 } else {
 ;the k ne kstart cases
 series_align,m2(i1:i2,j1:j2),ix,iy,5 }
 ;view for george data
 tv,ft(m2(i1:i2,j1:j2))
 ;2/9/93, the dxp, dyp are previous image shifts if any
 deltac(0,kmo)=2.0*[ix,iy]	;the 2.0 is for compressed version
 type,'k, offsets =',k, 2.0*ix, 2.0*iy
 
 ;some checks
 while (1) {
 sanity=3
 ;for the george data, we don't want shifts ! so set high limit
 if abs(ix) gt 10 or abs(iy) gt 10 then {
 sanity=sanity-1
 ;a large jump, if not too large we shift image and repeat CC
 printf,1,'k, LARGE offsets =',k,deltac(0,kmo), deltac(1,kmo)
	 if abs(ix) gt 40 or abs(iy) gt 40 or sanity lt 0 then {
	 ;really too big (actually 80 pixels because ix, iy in compressed space)
	 printf,1,'MAJOR PROBLEM, offset excceeds limit.'
	 ty,'MAJOR PROBLEM, offset excceeds limit.'
	 run align_reset
	 ;giving up on this image but hoping that the next one may still be
	 ;OK involves using the previous (good?) image for CC and
	 ;making delta=0, any value can be used for shiftc, make it 100 as
	 ;a flag, defer the untwist call if it was set
	 if #twist_option eq 0 then #twist_option=1
	 deltac(0,kmo)=[0.0,0.0]
	 shiftc(0,kmo)=[100.,100.]
	 break	}
 ;restore the previous transform into $a, $b, etc
 run align_reset
 ;try again with a shift, include old shift if any
 dx=rfix(deltac(0,kmo)+shiftc(0,kmo))	dy=rfix(deltac(1,kmo)+shiftc(1,kmo))
 shiftc(0,kmo)=[dx,dy]
 ii1=i1-dx	ii2=i2-dx	jj1=j1-dy	jj2=j2-dy
 ;need to check limits
 if ii1 lt 0 or jj1 lt 0 or ii2 ge nx or jj2 ge ny then {
	 printf,1,'MAJOR PROBLEM, subarea shift hits edge'
	 ty,'MAJOR PROBLEM, subarea shift hits edge'
	 ;same situation as above, hope for better images later on
	 if #twist_option eq 0 then {
	 	#twist_option=1
		type,'ATTENTION: #twist_option now set to 1'
		printf,1,'ATTENTION: #twist_option now set to 1'
		}
	 deltac(0,kmo)=[0.0,0.0]
	 shiftc(0,kmo)=[100.,100.]
	 break	}
 series_align,m2(ii1:ii2,jj1:jj2),ix,iy,5
 deltac(0,kmo)=2.0*[ix,iy]
 type,'new offsets =',deltac(0,kmo),deltac(1,kmo) -
	 ,' shifts =',shiftc(0,kmo),shiftc(1,kmo)
 
 } else break		;when ix and iy small enough, break out of while
 }
 
 }	;of the image loop
 ;now include the image shifts in the deltac's
 deltac=deltac+shiftc
 deltac(0,1)=deltac(*,1:*)-shiftc(*,0:(nf-2))
 endsubr
 ;==============================================================================
block align_reset
 switch,$a2,$a	switch,$b2,$b	switch,$c2,$c	switch,$d2,$d
 $align_count-=1
 endblock
 ;==============================================================================
subr dhipass(d,dxs,dys)
 ;conditions the offsets for jitter removal
 ;this used to be included in untwist
 ;uses results from getalign
 #max_shift=3
 dx=d(0,*)
 dy=d(1,*)
 dxs=runsum(dx)	dys=runsum(dy)
 ;drift options, #drift=0 for rigid, 1 for linear drift, gt 1 for keeping
 ;a smoothed version of the displacement with boxcar width = #drift
 if #drift eq 1 then { dxs=detrend(dxs,1)	dys=detrend(dys,1)  }
 ;2/9/93 gsmooth (gaussian smooth) used in place of boxcar
 ;use #drift as the fwhm of the gaussian
 if #drift gt 1 then {
	 dxs=ghipass(dxs, #drift, #max_shift)
	 dys=ghipass(dys, #drift, #max_shift)
 }
 ;11/2/90 - may be safer to use mean as reference, avoids single glitch problem
 dxs=dxs-mean(dxs)		dys=dys-mean(dys)
 endsubr
 ;==============================================================================
subr untwist(kstart,kend,name_in,name_out,r_ref,r,dxs,dys,s1,s2)
 ;2/3/93 subroutine version adapted from block version in alignmovies
 ;remove rotation and previously computed rigid alignment
 ;get the rigid alignments to apply to COR files
 
 if #twist_option ne 1 then {		;actually do it for options 0, 3, and 2
 
 ;-----------------loop over the images now ------------------------
 for k=kstart,kend do {
 kmo=k-kstart
 fseek,x,fns(name_in,k),head
 ndisk=#last_disk	;save for possible delete
 ;for the first value of k, check the array size
 if k eq kstart then { nx=dimen(x,0) ny=dimen(x,1) xrc=nx/2 yrc=ny/2 }
 m2=rotate3(x,r(kmo),xrc,yrc,dxs(kmo),dys(kmo))
 h2=head+smap(byte(13))+'r_ref ='+string(r_ref)
 if num_elem(h2) gt 256 then h2=h2(0:255)
 fcwrite,m2,fns(name_out,k),h2
 ty,'wrote: ',fns(name_out,k)
 
 ;get sample images for kstart and kend if #twist_option = 0 (case = 1 below)
 if #twist_option eq 0 then {
 if k eq kstart then s1=ft(compress(m2,4))
 if k eq kend then s2=ft(compress(m2,4))
 }
 
 ;the following deletes input file if #dflag = 1
 if #dflag eq 1 then {
 sdel='rm '+eval('#disk'+istring(ndisk,2,2))+fns(name_in,k)
 spawn,sdel	}	;end of delete option
 ;memstat		;ck memory status
 
 }	;of the image loop
 } else {				;for option 1, just make 2 samples
 
 k=kstart
 fseek,x,fns(name_in,k),head
 kmo=k-kstart
 nx=dimen(x,0) ny=dimen(x,1) xrc=nx/2 yrc=ny/2
 m2=rotate3(x,r(kmo),xrc,yrc,dxs(kmo),dys(kmo))
 ;get sample images for kstart and kend
 s1=ft(compress(m2,4))
 k=kend
 fseek,x,fns(name_in,k),head
 kmo=k-kstart
 m2=rotate3(x,r(kmo),xrc,yrc,dxs(kmo),dys(kmo))
 s2=ft(compress(m2,4))
 
 }
 endsubr
 ;=============================================================================
subr getstretch,kstart,kend,name_in,delta,grid_set,grid_clips
 ;get final destretch grid size
 nf=kend-kstart+1
 ngx=grid_set(dimen(grid_set)-1)		ngy=ngx
 ty,'ngx final =', ngx
 ;the array to hold the grid results, make it a zero scalar for the moment
 delta=0
 m1=0	m2=0	bad=0
 ;important, read only the best of each pair
 fzread,best,'/pore2/people/shine/george.best'
 
 for k=kstart,kend do {
 kmo=k-kstart
 t1=!systime
 ;special for run with George (Simon/Brandt) data
 
 kbest=best(k-1)
 ty,'image #',kbest
 fseek,xfull,fns(name_in,kbest),head
 x=xfull(10:709,10:709)		;for the 7 hour segment
 ;fseek,xfull,fns(name_in,kbest),head		;read in the image
 ;x=xfull(*,0:369)	;subarea of subarea, (1160:1849,1041:1510) in original
 
 nq=dimen(x,0)
 if #trim_edge eq 1 then {
 case
 nq eq 512:m2=x(16:495,16:495)	;just a little off the edges
 nq eq 1024:m2=x(32:991,32:991)	;the 1024 case
 else:switch,m2,x		;others, no cutting, handy for restretching
 endcase
 } else { switch,m2,x  }
 ;if this is not the first, do the alignment
 if k ne kstart then {
 dq=dsgridnest(m1,m2,grid_set,grid_clips)
 ;define the output array if we haven't yet
 if isarray(delta) eq 0 then {
  delta=zero(fltarr(2,dimen(dq,1),dimen(dq,2),nf))
 }
 ;zero any accumulated offsets > #max_stretch_shift
 in=sieve( abs(dq) gt #max_stretch_shift )
 if isarray(in) then { dq(in)=0.0   iq = num_elem(in)  !badmatch += iq
   ty,'extra (accumulated offset > limit) bad matches =',iq
  printf,1,'extra (accumulated offset > limit) bad matches =',iq
   }
 delta(0,0,0,kmo)=dq
 xq=mean(abs(dq(0,*,*)))	yq=mean(abs(dq(1,*,*)))		zq=max(abs(dq))
 t,'mean dx,dy =',xq,yq,' max =',zq, ' !badmatch =',!badmatch
 printf,1,'mean dx,dy =',xq,yq,' max =',zq, ' !badmatch =',!badmatch
 bad=bad+!badmatch
 ;note that m2 is unchanged (it will be the next m1)
 }	;of the k ne 0 conditional
 switch,m1,m2
 t2=!systime
 type,'the clock time for this image was',t2-t1,k
 }	;of the image loop
 printf,1,'!badmatch total =',bad
 endsubr
 ;=======================================================================
subr destretch,kstart,kend,name_in,name_out,delta
 !stretchmark=1
 ;IMPORTANT NOTE - nsm controls the unsharp filtering of the LCT results
 ;it effectively high pass filters the getstretch results to remove the
 ;short time scale seeing variations from the movies while retaining flows
 ;on longer time scales, good values for nsm are dependent on the data set
 nsm=#drift
 t1=!systime
 delx=delta(0,*,*,*)	dely=delta(1,*,*,*)	delta=0	;last to save memory
 ngx=dimen(delx,0)	ngy=dimen(delx,1)
 ;smooth the displacements and various other tricks
 ;put the time dimension in the first dimension
 ;to make the running sum and detrend easier
 delx=delx(>1,>2,>0)	dely=dely(>1,>2,>0)
 nsm=fix(nsm)
 nt=dimen(delx,0)
 ;loop over the grid, assumes no awful glitches, don't smooth if nsm=0
 ;2/19/93, linear trends are removed first
 for j=0,ngy-1 do for i=0,ngx-1 do {
 xq=runsum(delx(*,i,j))
 if nsm gt 0 then { delx(0,i,j) = ghipass(xq-trend(xq,1), nsm, 3)
	 } else { delx(0,i,j) = xq - mean(xq) }
 xq=runsum(dely(*,i,j))
 if nsm gt 0 then { dely(0,i,j) = ghipass(xq-trend(xq,1), nsm, 3)
	 } else { dely(0,i,j) = xq - mean(xq) }
 }
 ;now put the dimensions back in the original order
 delx=delx(>2,>0,>1)	dely=dely(>2,>0,>1)
 ;reconstruct the 4-D delta
 delta=[ [delx],[dely] ]
 delta=delta(>1,>2,>3,>0)
 m2=0
 t2=!systime
 ty,'clock time to process offsets =',t2-t1
 printf,1,'clock time to process offsets =',t2-t1
 ;important, read only the best of each pair
 fzread,best,'/pore2/people/shine/george.best'

 ;-----------------------------------------------------------------------------
 for k=kstart,kend do {
 kmo=k-kstart
 t1=!cputime
 ;fseek,x,fns(name_in,k),head		;read in the image
 ;special for run with George (Simon/Brandt) data
 
 kbest=best(k-1)
 ty,'image #',kbest
 fseek,xfull,fns(name_in,kbest),head
 x=xfull(10:709,10:709)		;for the 7 hour segment
 ;x=xfull(*,0:369) ;11 hour subarea of subarea, (1160:1849,1041:1510) in original
 ndisk=#last_disk	;save for possible delete
 nq=dimen(x,0)
 if #trim_edge eq 1 then {
 case
 nq eq 512:m2=x(16:495,16:495)	;just a little off the edges
 nq eq 1024:m2=x(32:991,32:991)	;the 1024 case
 else:switch,m2,x		;others, no cutting
 endcase
 } else { switch,m2,x }
 mq=stretch(m2,delta(*,*,*,kmo))
 fcwrite,mq,fns(name_out,k),head
 ;the following deletes input file if #dflag = 1
 if #dflag eq 1 then {
 sdel='rm '+eval('#disk'+istring(ndisk,2,2))+fns(name_in,k)
 spawn,sdel	}	;end of delete option
 t2=!cputime
 t,'the cpu time for this image was',t2-t1,k
 }	;of the image loop
 endsubr
 ;=======================================================================