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
 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)
 ;4/6/96 make compressing by factor of 2 an option, it was default, want
 ;to be able to have more accurate values even thought it takes longer
 ;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
 if defined(#ccx2flag) eq 0 then #ccx2flag=0	;no longer compress by default
 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 {
 t1 = !systime
 kmo=k-kstart
 fseek,x,fns(name_in,k),head		;read in the image
 ;get sample images for kstart and kend
 
 if k eq kstart then s1=ft(compress(x,4))
 if k eq kend then s2=ft(compress(x,4))
 
 ;for the first value of k, set up subarea parameters
 if k eq kstart then { 
   nx=dimen(x,0)  ny=dimen(x,1)
   xrc=nx/2	yrc=ny/2
   if #usersubarea eq 0 then {	;default choices if not user specified
   nq=nx<ny  nsub=384
   if nq gt 768 then {	nsub=768 }
   ;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
   xrc = (i1+i2)/2	yrc = (j1+j2)/2
   }	;end of #usersubarea conditional
   t,'selected subarea ranges are',i1,i2,j1,j2
   t,'rotation center =', xrc, yrc
   printf,1,'selected subarea ranges are',i1,i2,j1,j2	}
 ;2/17/93 - seeing measure added, not very efficient, consider using
 ;something from transform space since we do an fft anyhow
 
 xq=x(i1:i2,j1:j2)
 meanint(kmo)=mean(xq)			;mean intensity
 seeing(kmo)=mean(abs(differ(xq)))/meanint(kmo)	;"seeing" measure
 
 ;2/7/97  we should be rotating about the subarea, we can't just
 ;take subarea and then rotate (corners get cut) so we rotate entire
 ;image about the center of subarea, for small subareas this is inefficient
 ;rotate ?
 if #rotate ne 0 then m2=rotate(x,r(kmo),xrc,yrc) else switch,m2,x
 ;use a compressed (by 2) version of the data, 11/5/89
 ;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
 m2s = m2(i1:i2,j1:j2)
 if #ccx2flag eq 1 then m2c = compress(m2s,2) else switch, m2s,m2c
 start_align, m2c, lowcut, hicut
 ix=0	iy=0
 } else {
 ;the k ne kstart cases
 m2s = m2(i1:i2,j1:j2)
 if #ccx2flag eq 1 then m2c = compress(m2s,2) else switch, m2s,m2c
 series_align, m2c, ix, iy, 5
 }
 tv,ft(m2c),0,0,0
 ;2/9/93, the dxp, dyp are previous image shifts if any
 if #ccx2flag eq 1 then {
  deltac(0,kmo)=2.0*[ix,iy]	;the 2.0 is for compressed version
  type,'k, offsets =',k, 2.0*ix, 2.0*iy
 } else {
  deltac(0,kmo)=[ix,iy]
  type,'k, offsets =',k, ix, iy
 }
 ;some checks
 sanity=5
 while (1) {
 if abs(ix) gt 2 or abs(iy) gt 2 then {
 sanity=sanity-1
 ty,'sanity =', sanity
 ;a large jump, if not too large we shift image and repeat CC
 ty,'k, LARGE offsets =',k,deltac(0,kmo), deltac(1,kmo)
 printf,1,'k, LARGE offsets =',k,deltac(0,kmo), deltac(1,kmo)
   if abs(ix) gt 60 or abs(iy) gt 60 or sanity lt 0 then {
   ;really too big (actually 120 pixels if 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)=[1000.,1000.]
   break	}
 ;restore the previous transform into $a, $b, etc, and save current $a2, 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	}
 m2s = m2(ii1:ii2,jj1:jj2)
 ty,'shifted cutout =', ii1,ii2,jj1,jj2
 if #ccx2flag eq 1 then m2c = compress(m2s,2) else switch, m2s,m2c
 tv,ft(m2c),0,0,0
 series_align, m2c, ix, iy, 5
 if #ccx2flag eq 1 then {
  deltac(0,kmo)=2.0*[ix,iy]	;the 2.0 is for compressed version
  type,'k, offsets =',k, 2.0*ix, 2.0*iy
 } else {
  deltac(0,kmo)=[ix,iy]
  type,'k, offsets =',k, ix, iy
 }
 type,'new offsets =',deltac(0,kmo),deltac(1,kmo) -
	 ,' shifts =',shiftc(0,kmo),shiftc(1,kmo)
 ;9/24/94, however, to avoid propagating the shifts forever, we still need
 ;the result for the unshifted m2 in the hopes that the next one is
 ;aligned with it, so put the old $a2, etc saved as $a3 into the $a and friends
 switch,$a,$a3	switch,$b,$b3	switch,$c,$c3	switch,$d,$d3 
 $a3=0	$b3=0	$c3=0	$d3=0
 } else break		;when ix and iy small enough, break out of while
 }
 t2 = !systime
 ty,'time = ', t2-t1
 }	;of the image loop
 ;now include the image shifts in the deltac's
 deltac=deltac+shiftc
 ;don't think the next line is needed now that the shift has been removed
 ;in the CC
 ;;;deltac(0,1)=deltac(*,1:*)-shiftc(*,0:(nf-2))
 endsubr
 ;==============================================================================
block align_reset
 ;situation, we want to do the CC over, the new values are in $a and the
 ;old ones in $a2 (destined to be overwritten), we put the old values back
 ;in $a
 switch,$a2,$a	switch,$b2,$b	switch,$c2,$c	switch,$d2,$d
 $align_count-=1	;also knock down the count (or parity is wrong)
 ;we may also want to save those new values (now in $a2) so put them in a
 ;series $a3, etc
 switch,$a2,$a3	switch,$b2,$b3	switch,$c2,$c3	switch,$d2,$d3 
 endblock
 ;==============================================================================
subr getrigid_fl(flist,r,deltac,seeing,meanint,s1,s2)
 ;5/23/2000 - some restart situations added
 ;this fl version uses file list
 ;4/6/96 make compressing by factor of 2 an option, it was default, want
 ;to be able to have more accurate values even thought it takes longer
 ;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
 
 zeroifundefined, #twist_option, #usersubarea, #ccx2flag
 ;no longer compress by default
 nfiles = num_elem(flist)
 nf = nfiles
 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
 restart_flag = 0
 bad_series = 0
 ;loop over the images now
 
 ;-----------------------------------------------------------------------------
 kstart = 0
 kend = nfiles-1
 for k = kstart, kend do {
 t1 = !systime
 kmo=k-kstart
 fseek,x,flist(k),head		;read in the image
 ;get sample images for kstart and kend
 
 if k eq kstart then s1=ft(compress(x,4))
 if k eq kend then s2=ft(compress(x,4))
 
 ;for the first value of k, set up subarea parameters
 if k eq kstart or restart_flag then { 
   nx=dimen(x,0)  ny=dimen(x,1)
   xrc=nx/2	yrc=ny/2
   if #usersubarea eq 0 then {	;default choices if not user specified
   nq=nx<ny  nsub=384
   if nq gt 768 then {	nsub=768 }
   ;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
   xrc = (i1+i2)/2	yrc = (j1+j2)/2
   }	;end of #usersubarea conditional
   t,'selected subarea ranges are',i1,i2,j1,j2
   t,'rotation center =', xrc, yrc
   printf,1,'selected subarea ranges are',i1,i2,j1,j2
   }

 ;2/17/93 - seeing measure added, not very efficient, consider using
 ;something from transform space since we do an fft anyhow
 
 xq=x(i1:i2,j1:j2)
 meanint(kmo)=mean(xq)			;mean intensity
 seeing(kmo)=mean(abs(differ(xq)))/meanint(kmo)	;"seeing" measure
 
 ;2/7/97  we should be rotating about the subarea, we can't just
 ;take subarea and then rotate (corners get cut) so we rotate entire
 ;image about the center of subarea, for small subareas this is inefficient
 ;rotate ?
 if #rotate ne 0 then m2=rotate(x,r(kmo),xrc,yrc) else switch,m2,x
 ;use a compressed (by 2) version of the data, 11/5/89

 ;if this is not the first, do the alignment
 if k eq kstart or restart_flag then { 
 restart_flag = 0
 ;for the first one we need to set up alignment code
 m2s = m2(i1:i2,j1:j2)
 if #ccx2flag eq 1 then m2c = compress(m2s,2) else switch, m2s,m2c
 start_align, m2c, lowcut, hicut
 ix=0	iy=0
 } else {
 ;the k ne kstart cases
 m2s = m2(i1:i2,j1:j2)
 if #ccx2flag eq 1 then m2c = compress(m2s,2) else switch, m2s,m2c
 series_align, m2c, ix, iy, 5
 }
 tv,ft(m2c),0,0,0
 ;2/9/93, the dxp, dyp are previous image shifts if any
 if #ccx2flag eq 1 then {
  deltac(0,kmo)=2.0*[ix,iy]	;the 2.0 is for compressed version
  type,'k, offsets =',k, 2.0*ix, 2.0*iy
 } else {
  deltac(0,kmo)=[ix,iy]
  type,'k, offsets =',k, ix, iy
 }
 ;some checks
 sanity=5
 while (1) {
 if abs(ix) gt 2 or abs(iy) gt 2 then {
 sanity=sanity-1
 ty,'sanity =', sanity
 ;a large jump, if not too large we shift image and repeat CC
 ty,'k, LARGE offsets =',k,deltac(0,kmo), deltac(1,kmo)
 printf,1,'k, LARGE offsets =',k,deltac(0,kmo), deltac(1,kmo)
   if abs(ix) gt 200 or abs(iy) gt 200 or sanity lt 0 then {
   ;really too big
   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.]
   bad_series += 1
   ;if this happens twice in a row, probably a permanent shift
   if bad_series ge 2 then {
   	restart_flag = 1  bad_series = 0
   	ty,'trying a reset, expect some bad offsets'
   	}
   break	}
   bad_series = 0
 ;restore the previous transform into $a, $b, etc, and save current $a2, 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)=[1000.,1000.]
   break	}
 m2s = m2(ii1:ii2,jj1:jj2)
 ty,'shifted cutout =', ii1,ii2,jj1,jj2
 if #ccx2flag eq 1 then m2c = compress(m2s,2) else switch, m2s,m2c
 tv,ft(m2c),0,0,0
 series_align, m2c, ix, iy, 5
 if #ccx2flag eq 1 then {
  deltac(0,kmo)=2.0*[ix,iy]	;the 2.0 is for compressed version
  type,'k, offsets =',k, 2.0*ix, 2.0*iy
 } else {
  deltac(0,kmo)=[ix,iy]
  type,'k, offsets =',k, ix, iy
 }
 type,'new offsets =',deltac(0,kmo),deltac(1,kmo) -
	 ,' shifts =',shiftc(0,kmo),shiftc(1,kmo)
 ;9/24/94, however, to avoid propagating the shifts forever, we still need
 ;the result for the unshifted m2 in the hopes that the next one is
 ;aligned with it, so put the old $a2, etc saved as $a3 into the $a and friends
 switch,$a,$a3	switch,$b,$b3	switch,$c,$c3	switch,$d,$d3 
 $a3=0	$b3=0	$c3=0	$d3=0
 } else break		;when ix and iy small enough, break out of while
 }
 t2 = !systime
 ty,'time = ', t2-t1
 }	;of the image loop
 ;now include the image shifts in the deltac's
 deltac=deltac+shiftc
 ;don't think the next line is needed now that the shift has been removed
 ;in the CC
 ;;;deltac(0,1)=deltac(*,1:*)-shiftc(*,0:(nf-2))
 endsubr
 ;==============================================================================
subr dhipass(d,dxs,dys)
 ;conditions the offsets for jitter removal
 ;this used to be included in untwist
 ;uses results from getalign
 #max_shift=10
 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 gaussian fwhm = #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,kstartoffsets)
 ;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 - kstartoffsets
 fseek,x,fns(name_in,k),head
 ndisk=#last_disk	;save for possible delete
 ;if a w6302m60rcp, then patch
 ;if strpos(name_in,'w6302m60rcp') gt 1 then {  }
 ;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 !crunch_slice=5
 if symdtype(x) eq 0 then !crunch_slice=1
 }
 m2=rotate3(x,r(kmo),xrc,yrc,dxs(kmo),dys(kmo))
 
 ;if input was fits, let it be so for output
 if $last_type eq 1 then {
 wrfitscrun, m2, fns(name_out, k), head
 } else {
 h2=head+smap(byte(13))+'r_ref ='+string(r_ref)
 if num_elem(h2) gt 256 then h2=h2(0:255)
 if symdtype(m2) le 2 then { 	fcwrite,m2,fns(name_out,k),h2 } else {
 				fzwrite,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
 ty,'deleted'	}	;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 - kstartoffsets
 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- kstartoffsets
 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,seeing
 ;1/5/94 upgraded for rectangular arrays
 ;get final destretch grid size
 nf=kend-kstart+1
 ngx=grid_set(dimen(grid_set)-1)
 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
 !match_messages = 0		;suppress all the messages
 for k=kstart,kend do {
 kmo=k-kstart
 t1=!systime
 fseek,x,fns(name_in,k),head		;read in the image
 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
 nq eq 1280:m2=x(32:1247,32:991)	;the 1280 by 1024 case (probably)
 nq eq 1534:m2=x(32:1501,32:991)
 nq eq 1536:m2=x(32:1501,32:991)
 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,kstartgrid
 !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
 ;sometimes the conditioning takes a long time so put some indicators
 ;in every 10 minutes
 tlast=!systime
 delx=delx(>1,>2,>0)
 if (!systime - tlast) ge 600. then {ty,'done with x transpose' tlast=!systime}
 dely=dely(>1,>2,>0)
 if (!systime - tlast) ge 600. then {ty,'done with y transpose' tlast=!systime}
 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 {
 if (!systime - tlast) ge 600. then {ty,'at i,j=',i,j  tlast=!systime}
 xq=runsum(delx(*,i,j))
 if nsm gt 0 then { delx(0,i,j) = ghipass(xq-trend(xq,1), nsm, 8)
	 } 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, 8)
	 } else { dely(0,i,j) = xq - mean(xq) }
 }
 ;now put the dimensions back in the original order
 delx=delx(>2,>0,>1)
 if (!systime - tlast) ge 600. then {ty,'back x transpose' tlast=!systime}
 dely=dely(>2,>0,>1)
 if (!systime - tlast) ge 600. then {ty,'back y transpose' tlast=!systime}
 ;reconstruct the 4-D delta
 delta=[ [delx],[dely] ]
 if (!systime - tlast) ge 600. then {ty,'concat' tlast=!systime}
 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
 ;-----------------------------------------------------------------------------
 for k=kstart,kend do {
 kmo=k-kstartgrid
 t1=!cputime
 fseek,x,fns(name_in,k),head		;read in the image
 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
 nq eq 1280:m2=x(32:1247,32:991)	;the 1280 by 1024 case (probably)
 nq eq 1534:m2=x(32:1501,32:991)
 nq eq 1536:m2=x(32:1501,32:991)
 else:switch,m2,x		;others, no cutting
 endcase
 } else { switch,m2,x }
 mq=stretch(m2,delta(*,*,*,kmo))
 if symdtype(m2) le 2 then { 	fcwrite,mq,fns(name_out,k),head } else {
 				fzwrite,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
 ;=======================================================================