;for getting transforms of movies too big for memory
;this version is specialized for 200x200x128
;it could be (carefully!) modified for other applications
@READNEXT
;=========================================================================
block bigpass1
;currenly hardwired for pore region, 128 frames and a 200x200 subarea
;6/2/86  wired for sunspot series
name='DSK6:[soup]C20N038.RCC'	;first name
;name='DSK5:[FERGUSON]C00438.RCC'	;first name
;6/5/86 starting sunspot data with 438 to get interesting things at end
;6/6/86 also starting C20 series at 38 to get end
CLOSE,2
blkio,2,'SOUP:BIGPOWER.HOG'	;open our storage file
rw=assoc(2,fltarr(128))		;rw is one physical block in size
;note that for this pass, the file is divided into 4 parts, each of
;10,368 blocks to hold 4 128x101x101 floating point arrays
;during the second pass, this will be converted to 8 65x101x101 arrays
nxf=101
nyf=101	;appropiate for a 200x200 area
nxfm=nxf-1	nyfm=nyf-1
nt=128
ndt=32	;the number to block
ntout=nt/ndt	;the outer loop
nbt=1	;# of physical blocks for time index range
k=0	;time counter
for ktout=0,ntout-1 do begin
;reconstruct time buffers for each pass (they get transposed!)
b1=fltarr(nyf,nxf,ndt)
b2=fltarr(nyf,nxf,ndt)
b3=fltarr(nyf,nxf,ndt)
b4=fltarr(nyf,nxf,ndt)
for kt=0,ndt-1	do begin	;time buffer loop
type,'starting ',name
readnext,x,name
;image in x, get the 2-D transform (easy part)
;x=x(28:227,10:209)
x=x(28:227,28:227)
;6/1/86 - try a little temporal apodizing, a sine bell of quarter period 4
if k le 3 then x=x*sin(k*#pi/8.)
if k ge 124 then x=x*sin( (127-k)*#pi/8.)
xmean=mean(x)
type,'k= ',k,' mean =',xmean
x=x-xmean
sc,x,sx,cx
sx=sx(>1,>0)	;transpose for second pass
cx=cx(>1,>0)
sc,sx,sysx,cysx
sc,cx,sycx,cycx
;load the time buffers, note this also converts to single precision
b1(0,0,kt)=sysx
b2(0,0,kt)=cysx
b3(0,0,kt)=sycx
b4(0,0,kt)=cycx
nextname,name
k+=1	;bump k
end	;end of kt loop
;now transpose and unload buffers
b1=b1(>1,>2,>0)
b2=b2(>1,>2,>0)
b3=b3(>1,>2,>0)
b4=b4(>1,>2,>0)
iq=0
kq=ndt*ktout
type,'begin file output for ktout =',ktout
for i=0,nxfm do for j=0,nyfm do begin
;read in the designated block
buf=rw(iq)
buf(kq)=b1(*,j,i)	;load this part, don't change the rest
rw(iq)=buf		;write back to file
;that was buffer 1 (b1), now the same for the rest
buf=rw(iq+10368)
buf(kq)=b2(*,j,i)	;load this part, don't change the rest
rw(iq+10368)=buf		;write back to file
buf=rw(iq+20736)
buf(kq)=b3(*,j,i)	;load this part, don't change the rest
rw(iq+20736)=buf		;write back to file
buf=rw(iq+31104)
buf(kq)=b4(*,j,i)	;load this part, don't change the rest
rw(iq+31104)=buf		;write back to file
iq+=1	;bump iq
end	;of file write loop
end	;of ktout loop
close,1
X=0	SX=0	CX=0	SYSX=0	SYCX=0 CYSX=0 CYCX=0
endblock
;=========================================================================
block bigpass2
close,2
blkio,2,'SOUP:BIGPOWER.HOG'	;open our storage file
rw=assoc(2,fltarr(128))	;rw is again 1 block
;pass 1 left 4 large zones on the file, each is organized as a 128x101x101
;single precision array and separated by 10368 blocks
;these are already transposed for the temporal FFT, we only need to
;read and do a 1-D fft and store results
;
;zero some old buffers (if left over from pass 1) and set up new
b1=0	b2=0	b3=0	b4=0
b1=fltarr(663552)	;a bit larger than 65x101x101x4
b2=fltarr(663552)	;to fit 1/8 of file
rout=assoc(2,b1)
;
for k=0,3 do begin	;loop over the big sections
ty,'starting section ',k
iq=k*10368	;section bias
jq=0
for j=0,10200 do begin	;divide into 10201 pieces
buf=doub(rw(iq+j))
sc,buf,st,ct
;now load these into b1 and b2
b1(jq)=float(st)
b2(jq)=float(ct)
jq=jq+65
end	;of j loop
;output the 2 buffers
ty,'starting output of buffers'
rout(2*k)=b1
rout(2*k+1)=b2
end	;of k loop
;that should be it!
;results are stored in 8 sections, each is 65x101x101
;they are on 5184 block boundaries and contain 2652260 bytes each
;(a bit over 5180 blocks)
endblock
;=========================================================================
BLOCK BACKPASS1
;the reverse of BIGPASS2
;5/27/86  contains a built in wedge filter for test
;construct some parts of the filter, we try to save memory
;;;VT=26.5*INDGEN(FLTARR(65))
VT=18.24*INDGEN(FLTARR(65))	;for c20 data
DELS=1				;width of gaussian apodizer for velocity
IN=INDGEN(FLTARR(101*101))
IX=IN/101   IY=IN%101
IN=SQRT(IX*IX+IY*IY)>1
VS=1.0/IN
;a high pass filter in space only
;FSPACE=(IN GE 20.0)+(IN LT 20.)*EXP(-((IN-20)/10)^2)
;note that VS is the inverse of the symmetrized spatial index and can
;also be used for spatial filter cutoffs
IX=0	IY=0	IN=0
close,2
blkio,2,'SOUP:BIGPOWER.HOG'	;open our storage file
rout=assoc(2,fltarr(663552)
rw=assoc(2,fltarr(128))	;rw is again 1 block
for k=0,3 do begin	;loop over the big sections
ty,'starting section ',k
;input the 2 buffers
b1=0	b2=0	;first free unused area
b1=rout(2*k)
b2=rout(2*k+1)
iq=k*10368	;section bias
jq=0
for j=0,10200 do begin	;divide into 10201 pieces
jqq=jq+64
;calculate the short filter for this pass
V=VT*VS(J)	;actual velocity, values are in km/s
;;;F=(V GE 6.0) AND (V LE 9.0)	;acoustic wedge
;;;F=(V GE 5.0) AND (V LE 8.5)	;acoustic wedge
;next 2 lines are the new apodized filter
;F=(V GT 8.0)*EXP(-((V-8.0)^2)/DELS)+(V LT 5.5)*EXP(-((V-5.5)^2)/DELS)
;F=F+( (V LE 8.0) AND (V GE 5.5) )
;;;F=V LE 3.0	;subsonic cone
f=v LE 5.0	;larger subsonic cone
;;;F=V GT 10.0	;supersonic cone
;test without a filter
;;;f=zero(v)+1.0
;F=F*FSPACE(J)
st=doub(b1(jq:jqq))*F
ct=doub(b2(jq:jqq))*F	;that should apply the wedge filter
scb,buf,st,ct
;now load into file
rw(iq+j)=float(buf)
jq=jq+65
end	;of j loop
end	;of k loop
close,2
;that should be it!
;results are stored in 4 sections, each is 128x101x101
;should be the same situation as at the end of forward pass 1
FSPACE=0	V=0	F=0	VT=0	VS=0
endblock
;=========================================================================
block backpass2
;reverse of bigpass1, we
;leave result in DSK5:[shine]result.3d
close,1
;;openu,1,'dsk5:[shine]result.3d'
openu,1,'dsk6:[soup]result.3d'
w=assoc(1,intarr(200,200))
CLOSE,2
blkio,2,'SOUP:BIGPOWER.HOG'	;open our storage file
rw=assoc(2,fltarr(128))		;rw is one physical block in size
nxf=101
nyf=101	;appropiate for a 200x200 area
nxfm=nxf-1	nyfm=nyf-1
nt=128
ndt=32	;the number to block
ntout=nt/ndt	;the outer loop
nbt=1	;# of physical blocks for time index range
k=0	;time counter
b1=0	b2=0	b3=0	b4=0
for ktout=0,ntout-1 do begin
;reconstruct time buffers for each pass (because we transpose them)
b1=fltarr(ndt,nyf,nxf)
b2=fltarr(ndt,nyf,nxf)
b3=fltarr(ndt,nyf,nxf)
b4=fltarr(ndt,nyf,nxf)
iq=0
kq=ndt*ktout
kq2=kq+ndt-1
type,'begin file input for ktout =',ktout
for i=0,nxfm do for j=0,nyfm do begin
;read in the designated block
buf=rw(iq)
b1(0,j,i)=buf(kq:kq2)
;that was buffer 1 (b1), now the same for the rest
buf=rw(iq+10368)
b2(0,j,i)=buf(kq:kq2)
buf=rw(iq+20736)
b3(0,j,i)=buf(kq:kq2)
buf=rw(iq+31104)
b4(0,j,i)=buf(kq:kq2)
iq+=1	;bump iq
end	;of file read loop
;now transpose (reverse of bigpass1)
b1=b1(>2,>0,>1)
b2=b2(>2,>0,>1)
b3=b3(>2,>0,>1)
b4=b4(>2,>0,>1)
for kt=0,ndt-1	do begin	;time buffer loop
;regenerate the sin/cosine arrays
sysx=b1(*,*,kt)
cysx=b2(*,*,kt)
sycx=b3(*,*,kt)
cycx=b4(*,*,kt)
scb,sx,sysx,cysx
scb,cx,sycx,cycx
sx=sx(>1,>0)	;transpose for second pass
cx=cx(>1,>0)
scb,x,sx,cx
xmin=min(x)	xmax=max(x)
type,'k =',k,' min, max =',xmin,xmax
w(k)=word(x)
k+=1	;bump k
end	;end of kt loop
end	;of ktout loop
close,1
close,2
endblock
;=========================================================================
BLOCK BIGPOWER
;assumes that bigpass1 and bigpass2 have been run and result is in
;BIGPOWER.HOG
close,2
blkio,2,'SOUP:BIGPOWER.HOG'	;open our storage file
rout=assoc(2,fltarr(663552)
sum=zero(fltarr(663065))
for k=0,7 do begin	;loop over the 8 pieces
b1=rout(k)
b2=b1(0:663064)
b1=0
sum=sum+b2*b2
b2=0
end
sum3=fltarr(65,101,101)
sum3(0)=sum
;the sum of the squares is now in sum3, this is 3-D power spectrum
;with direction summed
;
;now symmetrize the spatial dimensions
com3to2,sum3,p2
;p2 contains the 2D result, store it in DSK5:[SHINE]POWER.3D2
store,p2,'DSK5:[SHINE]POWER.3D2'
close,2
ENDBLOCK