;for getting transforms of movies using large amount of virtual memory
;this version is specialized for 270x270x100
;and assumes input of 256x256x83
;it was modified from FFT3DMEM.ANA
;it could be (carefully!) modified for other applications
;=========================================================================
subr nextname,name
;update the name
;2/18/87 - slightly upgraded so that it finds last rather than first period
;in name, thus avoid sub-directory defs.
;note that this is a very simple name updating and only affects the last 3
;chars. in a file name; e.g., input C20N003.RCC and get C20N004.RCC
temp=smap(reverse(bmap(name
iperiod=num_elem(temp)-strpos(temp,'.')-1
if iperiod ge 4 then {
;found the period, assume standard file numbering convention
iq=long(name(iperiod-3:iperiod-1))+1	;bump to next one
name=name(0:iperiod-4)+istring(iq,3,1)+name(iperiod:*)
}
endsubr
;=========================================================================
subr clean,x1,x2,x3,factor
;attempts to clean up specks on x2 by checking against x1
;dust candidates  are replaced with the mean of x1 and x3 at that position
mdust1=x2 gt factor*x1
mdust2=x2 lt x1/factor
ind=indgen(mdust1)
ic=sieve(ind,mdust1)
ic1=sieve(ic,x2(ic) gt factor*x3(ic))
ic=sieve(ind,mdust2)
ic2=sieve(ic,x2(ic) lt x3(ic)/factor)
if (num_dim(ic1) eq 0 and num_dim(ic2) eq 0) then begin
type,'no specks found'
end else begin
ic=[ic1,ic2]
nc=num_elem(ic)
type,'number of points =',nc
for k=0,nc-1 do begin
i=ic(k)
x2(i)=0.5*(x1(i)+x3(i))
end
end
endsubr
;=========================================================================
subr fftopen,name,name2
;open a pair of the eight trig arrays with names "name" and "name2"
;direct access with a record size of 51*4 bytes
openU,3,name,51*4
openU,4,name2,51*4
$w=assoc(3,fltarr(51,136))
$w2=assoc(4,fltarr(51,136))
endsubr
;=========================================================================
subr fftread,name,name2
;open a pair of the eight trig arrays with names "name" and "name2"
;direct access with a record size of 51*4 bytes
openr,3,name
openr,4,name2
$r=assoc(3,fltarr(51,136))
$r2=assoc(4,fltarr(51,136))
endsubr
;=========================================================================
subr fcread,x,name
;read and compress the next c47 image
f0read,x,name,hdr
;x=compress(x,2)	;only for C47Q not AQ or BQ
endsubr
;=========================================================================
block setupreads
;prime the system
;type,'starting ',name
fcread,x2,name
nextname,name
;type,'starting ',name
fcread,x3,name
nextname,name
;type,'starting ',name
fcread,x1,name
nextname,name
ifile=0
endblock
;=========================================================================
block getnext
ty,'ifile =',ifile
CASE
IFILE EQ 0:{	;already have 3
clean,x1,x2,x3,factor }
IFILE EQ 1:{ clean,x1,x2,x3,factor }
IFILE EQ nfil-1:{ clean,x1,x2,x3,factor }
ELSE:{
fcread,x3,name
nextname,name
clean,x1,x2,x3,factor }
ENDCASE
x16=float(x2)
;now redefine some arrays
SWITCH,X1,X2
SWITCH,X2,X3
ifile=ifile+1
endblock
;=========================================================================
block bigpass1
;name='QSA0:[soup]C47Q001.RCC'	;first name
name='QSA0:[soup]C47BQ001.RCC'	;first name
;setup the large buffers for the first pass
A=FLTARR(136,136,100)
B=FLTARR(136,136,100)
C=FLTARR(136,136,100)
D=FLTARR(136,136,100)
X=DBLARR(270,270)	;load with apodized image
AP=SIN(INDGEN(FLTARR(7))*#PI/14.)	;spatial skirt
APR=REVERSE(AP)
nfil=83
factor=1.2
RUN,SETUPREADS
RUN,GETNEXT
NT=100
klast=(100+nfil)/2 -1
ty,'klast =',klast
for k=0,nt-1 do begin
if (k gt 8 and k le klast) then { run,getnext }
type,'starting ',name,' k =',k
;image in X16
xmean=mean(X16)
type,' mean =',xmean
;pack X16 into larger FP X with window
for j=0,255 do X(7,j+7)=X16(*,j)
X=X-XMEAN
;now put a cosine bell skirt on it (note mean is now 0)
for j=0,6 do X(0,j)=X(*,7)*AP(j)
for j=263,269 do X(0,j)=X(*,262)*AP(269-j)
for j=0,269 do { X(0,j)=AP*X(7,j)    X(263,j)=APR*X(262,j)  }
;try a little temporal apodizing, a sine bell of quarter period 8
CASE
k lt 8:{  fac=sin(k*#pi/16) 	sc,fac*x,sx,cx }
k gt (klast+1):{  fac=sin( (99-k)*#pi/16)	sc,fac*x,sx,cx }
else { sc,x,sx,cx }
ENDCASE
sx=sx(>1,>0)	;transpose for second pass
cx=cx(>1,>0)
sc,sx,sysx,cysx
sc,cx,sycx,cycx
;load the coefficient buffers, note this also converts to single precision
D(0,0,k)=sysx
B(0,0,k)=cysx
C(0,0,k)=sycx
A(0,0,k)=cycx
end	;of k loop
;restore some memory
X=0	SX=0	CX=0	SYSX=0	SYCX=0 CYSX=0 CYCX=0
endblock
;=========================================================================
block midpass
;in this pass we do the final FFT along the time dimension
;then we apply a filter to it and immediately unFFT in time
;the spatial back FFT is in the next block
DELS=1				;width of gaussian apodizer for velocity
IN=INDGEN(FLTARR(136,136))
IX=IN/136   IY=IN%136
IN=SQRT(IX*IX+IY*IY)>1
VS=1.0/IN
;VT=61.6*indgen(fltarr(51))	;for C47 data, with 0.89 factor included
VT=30.8*indgen(fltarr(51))	;for C47AQ or BQ data, with 0.89 factor
VONE=ZERO(VT)+1.0
V=FLTARR(51,136)		;will be loaded
FS=FLTARR(51,136)
FDEX=INDGEN(LONG(FS))
;note that VS and VT are the space and time factors of velocity
;the 4 coeff. buffers are A,B,C, and D
;we identify each with X in turn in the following loop
x=0
for kb=0,3 do begin	;loop over the big sections
ncase kb
{ switch,x,a }
{ switch,x,a   switch,x,b  }
{ switch,x,b   switch,x,c  }
{ switch,x,c   switch,x,d  }	;got that all figured ?
endcase
for j=0,135 do begin
for i=0,135 do V(0,I)=VT*VS(I,J)	;velocity for this set
for i=0,135 do FS(0,I)=VONE/VS(I,J)	;an unclever way, depends only on I
;remember that logical results are either 0 or 1
F=V LE 5.0	;a subsonic filter
;filter constructed, do the final t transform (forward)
buf=X(*,j,*)
buf=buf(>1,>0)
if !fftdp eq 1 then buf=doub(buf)
sc,buf,st,ct
;now apply filter
ST=ST*F		CT=CT*F
scb,buf,st,ct	;back in time dimension
buf=float(buf)
buf=buf(>1,>0)
;putting it back a bit more difficult
for k=0,99 do X(0,j,k)=buf(*,k)
end	;of j loop
end	;of kb loop
switch,x,d	;restore d
;that should be it!
endblock
;=========================================================================
block fundmidpass
;specialized for a fundamental filter
;in this pass we do the final FFT along the time dimension
;then we apply a filter to it and immediately unFFT in time
;the spatial back FFT is in the next block
DELS=1				;width of gaussian apodizer for velocity
IN=INDGEN(FLTARR(136,136))
IX=IN/136   IY=IN%136
IN=SQRT(IX*IX+IY*IY)>1
VS=1.0/IN
;VT=61.6*indgen(fltarr(51))	;for C47Q data
VT=30.8*indgen(fltarr(51))	;for C47AQ or BQ data, with 0.89 factor
V=FLTARR(51,136)		;will be loaded
VTG=indgen(fltarr(51))
VTG=VTG*VTG	;square it
;VTG=0.1935*VTG		;for C47Q data (includes 0.89 factor)
VTG=0.0968*VTG		;for C47AQ or BQ data (includes 0.89 factor)
VONE=ZERO(VTG)+1.0
VG=FLTARR(51,136)		;will be loaded
FS=FLTARR(51,136)
FDEX=INDGEN(LONG(FS))
;g is 0.274 km/s/s
GRAVTOP=1.3*0.274
;note that VS and VTG are the space and time factors of acceleration (g)
;the 4 coeff. buffers are A,B,C, and D
;we identify each with X in turn in the following loop
x=0
for kb=0,3 do begin	;loop over the big sections
ncase kb
{ switch,x,a }
{ switch,x,a   switch,x,b  }
{ switch,x,b   switch,x,c  }
{ switch,x,c   switch,x,d  }	;got that all figured ?
endcase
for j=0,135 do begin
for i=0,135 do VG(0,I)=VTG*VS(I,J)	;velocity for this set
for i=0,135 do V(0,I)=VT*VS(I,J)	;velocity for this set
for i=0,135 do FS(0,I)=VONE/VS(I,J)	;an unclever way, depends only on I
;remember that logical results are either 0 or 1
;wedge with low cut and high cut
;F=(VG GE GRAVBOT) AND (VG LE GRAVTOP)	;fundamental movie
;everything above fundamental*1.3
;F=VG GT GRAVTOP 			;superfundamental
;everything below fundamental*0.7
F=VG LT GRAVBOT  			;subfundamental
;F=VG LT GRAVTOP  			;sub plus fundamental
;everything below fundamental*0.8 but above 4km/s
;;F=(VG LT GRAVBOT) AND (V GE 4.0)	;movie 113
;filter constructed, do the final t transform (forward)
buf=X(*,j,*)
buf=buf(>1,>0)
if !fftdp eq 1 then buf=doub(buf)
sc,buf,st,ct
;now apply filter
ST=ST*F		CT=CT*F
scb,buf,st,ct	;back in time dimension
buf=float(buf)
buf=buf(>1,>0)
;putting it back a bit more difficult
for k=0,99 do X(0,j,k)=buf(*,k)
end	;of j loop
end	;of kb loop
switch,x,d	;restore d
;that should be it!
endblock
;=========================================================================
BLOCK BACKPASS1
;the reverse of BIGPASS1
openu,1,'dsk6:[soup]result.3d'
w=assoc(1,intarr(270,270)
koffset=0
;koffset=93	;only for the BQ series, store in high half of result.3d
for k=0,nt-1 do begin
;unload the coefficient buffers, note this also converts to single precision
sysx=D(*,*,k)
cysx=B(*,*,k)
sycx=C(*,*,k)
cycx=A(*,*,k)
SCB,sx,sysx,cysx
SCB,cx,sycx,cycx
sx=sx(>1,>0)	;transpose for second pass
cx=cx(>1,>0)
SCB,x,sx,cx 
x=word(x)	;back to 16 bits
w(k+koffset)=x		;write out to file
end	;of k loop
close,1
endblock
;=========================================================================
BLOCK BIGPOWER
;assume that bigpass1 has been run (but not midpass!)
;the 4 coeff. buffers are A,B,C, and D
;compute the temporal FFT and sum the squares of all of the 8 results
SUM=ZERO(FLTARR(51,136,136))	;this will hold result
;we identify each with X in turn in the following loop
;the 4 arrays are not destroyed
x=0
for kb=0,3 do begin	;loop over the big sections
ncase kb
{ switch,x,a }
{ switch,x,a   switch,x,b  }
{ switch,x,b   switch,x,c  }
{ switch,x,c   switch,x,d  }	;got that all figured ?
endcase
for j=0,135 do begin
buf=X(*,j,*)
buf=buf(>1,>0)
if !fftdp eq 1 then buf=doub(buf)
sc,buf,st,ct
sum(0,0,j)=st*st+ct*ct+sum(*,*,j)
end
end	;of kb loop
switch,x,d	;restore d
;the sum of the squares is now in sum, this is 3-D power spectrum
;with direction summed
;
;now symmetrize the spatial dimensions
com3to2,sum,p2
;p2 contains the 2D result, store it in DSK5:[SHINE]POWER20CAL.3D2
store,p2,'DSK5:[SHINE]POWER20CAL.3D2'
close,2
ENDBLOCK
;=========================================================================
BLOCK BIG3DPOWER
;save entire 3-D power spectrum in 8 arrays and store on disk
;assume that bigpass1 has been run (but not midpass!)
;the 4 coeff. buffers are A,B,C, and D
;we identify each with X in turn in the following loop
;note that the arrays are destroyed
x=0
for kb=0,3 do begin	;loop over the big sections
ncase kb
{ switch,x,a }
{ switch,x,b  }
{ switch,x,c  }
{ switch,x,d  }	;got that all figured ?
endcase
ncase kb
{ fftopen,'dsk6:[soup]c20fft.ctcycx','dsk6:[soup]c20fft.stcycx' }
{ fftopen,'dsk6:[soup]c20fft.ctcysx','dsk6:[soup]c20fft.stcysx' }
{ fftopen,'dsk6:[soup]c20fft.ctsycx','dsk6:[soup]c20fft.stsycx' }
{ fftopen,'dsk6:[soup]c20fft.ctsysx','dsk6:[soup]c20fft.stsysx' }
endcase
for j=0,135 do begin
buf=X(*,j,*)
buf=buf(>1,>0)
if !fftdp eq 1 then buf=doub(buf)
sc,buf,st,ct
$w(j)=float(ct
$w2(j)=float(st
end
close,3		close,4
end	;of kb loop
ENDBLOCK
;============================================================================
BLOCK BACKFROMFFT
;assumes FFT is done and stored as 8 trig. arrays, each array is (51,136,136)
;this block does the reverse FFT in time and restores the 4 trig. arrays
;A,B,C,D  -- these are then used by BACKPASS1 to restore the data
;filters should be applied in this block before the reverse FFT
;
;the 4 coeff. buffers are A,B,C, and D
;we identify each with X in turn in the following loop
NT=100
A=0  B=0  C=0  D=0
;filter related stuff
IN=INDGEN(FLTARR(136,136))
IX=IN/136   IY=IN%136
IN=SQRT(IX*IX+IY*IY)>1
VS=1.0/IN
VT=16.41*indgen(fltarr(51))	;for C20 data
VONE=ZERO(VT)+1.0
V=FLTARR(51,136)		;will be loaded
FS=FLTARR(51,136)
FDEX=INDGEN(LONG(FS))
;note that VS and VT are the space and time factors of velocity
;
for kb=0,3 do begin	;loop over the big sections
x=fltarr(136,136,100)
ncase kb
{ fftread,'QSA0:[soup]c20fft.ctcycx','QSA0:[soup]c20fft.stcycx' }
{ fftread,'QSA0:[soup]c20fft.ctcysx','QSA0:[soup]c20fft.stcycx' }
{ fftread,'QSA0:[soup]c20fft.ctsycx','QSA0:[soup]c20fft.stsycx' }
{ fftread,'QSA0:[soup]c20fft.ctsysx','QSA0:[soup]c20fft.stsysx' }
endcase
for j=0,135 do begin
;filter related
for i=0,135 do V(0,I)=VT*VS(I,J)	;velocity for this set
;need the following line only for spatial boundaries
;for i=0,135 do FS(0,I)=VONE/VS(I,J)	;an unclever way, depends only on I
;construct the filter
F=V LE 4.0	;a subsonic filter
;end of filter related area
ct=$r(j)	;ct will be (51,136) array
st=$r2(j)
;apply filter
ct=f*ct		st=f*st
if !fftdp eq 1 then { ct=doub(ct)  st=doub(st) }
scb,buf,st,ct	;back in time dimension
buf=float(buf)
buf=buf(>1,>0)
;putting it back a bit more difficult
for k=0,99 do X(0,j,k)=buf(*,k)
end	;of j loop
close,3		close,4
ncase kb
{ switch,x,a }	;put x into a, etc, depending on kb
{ switch,x,b  }
{ switch,x,c  }
{ switch,x,d  }	;got that all figured ?
endcase
end	;of kb loop
ENDBLOCK