func trigmag(x,n)
 ;magnify an array using sin/cos fit
 sc2d,x,a,b,c,d
 ;now create transforms for a bigger image
 nxf = dimen(a,0)	nyf = dimen(a,1)
 nxf = (nxf-1)*n + 1	nyf = (nyf-1)*n + 1
 aa = zero(fltarr(nxf, nyf))
 bb = aa	cc = aa		dd = aa
 ;insert the tranforms of the original into these zeroed arrays
 insert, aa, a
 insert, bb, b
 insert, cc, c
 insert, dd, d
 sc2db,y,aa,bb,cc,dd
 return, y
 endfunc
 ;=========================================================================
func dctmag(x,n)
 ;magnify an array using a dct fit
 y= x
 dcqb, y
 y=y(>1,>0)
 dcqb, y
 ;now create transforms for a bigger image
 nxf = dimen(y,0)	nyf = dimen(y,1)
 y = y*(1./(nxf*8*nyf*8)
 aa = zero(fltarr(nxf*n, nyf*n))
 ;insert the tranforms of the original into these zeroed arrays
 insert, aa, y
 dcqf, aa
 aa=aa(>1,>0)
 dcqf, aa
 return, aa
 endfunc
 ;=========================================================================
subr sc2d,x,a,b,c,d
 ;given a 2-D array x, returns the sin/cosine coeffs. in a,b,c,d
 sc,x,sx,cx
 sx=sx(>1,>0)	;transpose for second pass
 cx=cx(>1,>0)
 sc,sx,b,c
 sx=0		;save memory
 sc,cx,d,a
 cx=0		;save memory
 endsubr
 ;=========================================================================
subr sc2db,x,a,b,c,d
 ;given the sin/cosine coeffs. in a,b,c, and d, returns a 2-D array x
 scb,sx,b,c
 scb,cx,d,a
 sx=sx(>1,>0)	;transpose for second pass
 cx=cx(>1,>0)
 scb,x,sx,cx
 sx=0		;save memory
 cx=0		;save memory
 endsubr
 ;=========================================================================
func power2d,x
 ;function version of power2d
 !fftdp=0
 sc2d,x,a,b,c,d
 xq=(a*a+b*b+c*c+d*d)
 a=2.*(c*d-a*b)
 b=0	c=0
 d=xq+a
 a=xq-a
 ;do the edges right
 ;since various terms were zero, this amounts to multiplying along the edges
 a(0,0)=4.*a(*,0)	;note kt is carried through
 n=dimen(a,0)-1
 m=dimen(a,1)-1
 a(0,m)=4.*a(*,m)
 for i=0,m do { a(0,i)=4.*a(0,i)   a(n,i)=4.*a(n,i)  }
 ;note that the corners are multiplied by 16
 ;now the same for d (although one edge overlaps a)
 d(0,0)=4.*d(*,0)	;note kt is carried through
 d(0,m)=4.*d(*,m)
 for i=0,m do { d(0,i)=4.*d(0,i)   d(n,i)=4.*d(n,i)  }
 ;now d is the upper rh block and a is the lower
 xq=0	;save memory
 ;construct P using reverses and concats
 ;also note that x and y indices are transposed in the a,b.. arrays
 ;hence we transpose P to get back the normal sense
 P=REVERSE(A,1)
 P=CONCAT(P,D(*,1:*))
 P=P(>1,>0)
 P=CONCAT(REVERSE(P),P(*,1:*))
 return,p
 endfunc
 ;=========================================================================
func power2dw, x, w
 ;a specialized version that leaves the sin/cos arrays in globals
 ;function version of power2d
 !fftdp=0
 sw = string(w)
 sc2d,x,a,b,c,d
 equate, eval('$data_ffta' + sw), a
 equate, eval('$data_fftb' + sw), b
 equate, eval('$data_fftc' + sw), c
 equate, eval('$data_fftd' + sw), d
 xq=(a*a+b*b+c*c+d*d)
 a=2.*(c*d-a*b)
 b=0	c=0
 d=xq+a
 a=xq-a
 ;do the edges right
 ;since various terms were zero, this amounts to multiplying along the edges
 a(0,0)=4.*a(*,0)	;note kt is carried through
 n=dimen(a,0)-1
 m=dimen(a,1)-1
 a(0,m)=4.*a(*,m)
 for i=0,m do { a(0,i)=4.*a(0,i)   a(n,i)=4.*a(n,i)  }
 ;note that the corners are multiplied by 16
 ;now the same for d (although one edge overlaps a)
 d(0,0)=4.*d(*,0)	;note kt is carried through
 d(0,m)=4.*d(*,m)
 for i=0,m do { d(0,i)=4.*d(0,i)   d(n,i)=4.*d(n,i)  }
 ;now d is the upper rh block and a is the lower
 xq=0	;save memory
 ;construct P using reverses and concats
 ;also note that x and y indices are transposed in the a,b.. arrays
 ;hence we transpose P to get back the normal sense
 P=REVERSE(A,1)
 P=CONCAT(P,D(*,1:*))
 P=P(>1,>0)
 P=CONCAT(REVERSE(P),P(*,1:*))
 a = 0		d = 0
 return,p
 endfunc
 ;=========================================================================
subr power2d,x,p
 !fftdp=0
 sc2d,x,a,b,c,d
 xq=(a*a+b*b+c*c+d*d)
 a=2.*(c*d-a*b)
 b=0	c=0
 d=xq+a
 a=xq-a
 ;do the edges right
 ;since various terms were zero, this amounts to multiplying along the edges
 a(0,0)=4.*a(*,0)	;note kt is carried through
 n=dimen(a,0)-1
 m=dimen(a,1)-1
 a(0,m)=4.*a(*,m)
 for i=0,m do { a(0,i)=4.*a(0,i)   a(n,i)=4.*a(n,i)  }
 ;note that the corners are multiplied by 16
 ;now the same for d (although one edge overlaps a)
 d(0,0)=4.*d(*,0)	;note kt is carried through
 d(0,m)=4.*d(*,m)
 for i=0,m do { d(0,i)=4.*d(0,i)   d(n,i)=4.*d(n,i)  }
 ;now d is the upper rh block and a is the lower
 xq=0	;save memory
 ;construct P using reverses and concats
 ;also note that x and y indices are transposed in the a,b.. arrays
 ;hence we transpose P to get back the normal sense
 P=REVERSE(A,1)
 P=CONCAT(P,D(*,1:*))
 P=P(>1,>0)
 P=CONCAT(REVERSE(P),P(*,1:*))
 end
 ;=========================================================================
subr dab,a,b,c,d,ix,iy
 ;puts up the 4 2D arrays
 ;and types means and maxes
 ty,'maxes for a,b,c, and d =',max(a),max(b),max(c),max(d)
 ;ty,'means for the same     =',mean(a),mean(b),mean(c),mean(d)
 nx=dimen(a,0) ny=dimen(a,1)
 ;load into larger array and strip out the zero columns and rows
 nd=2*nx-1	;note that nx should have been odd
 xq=fltarr(nd,nd)
 for j=0,ny-1 do xq(0,j)=a(*,j)
 for j=ny,nd-1 do xq(0,j)=c(*,j+1-ny)
 for j=0,ny-1 do xq(nx,j)=d(1:nx-2,j)
 for j=ny,nd-1 do xq(nx,j)=b(1:nx-1,j+1-ny)
 ty,'max, min of xq =',max(xq),min(xq)
 tv,xq,ix,iy
 endsubr
 ;=========================================================================
subr maskprep,x,nx,ny
 ;returns an array containing the freq. # for each element
 x=indgen(lonarr(nx,ny))
 y=x/nx
 x=x%nx
 x=sqrt(x*x+y*y)
 endsubr
 ;=========================================================================
func hipass(xin,k)
 ;returns a high pass filtered version of xin
 ;the cutoff freq. # is k
 sc2d,xin,a,b,c,d
 nx=dimen(a,0)
 ny=dimen(a,1)
 maskprep,xm,nx,ny	;get array of freq. #'s
 xm=xm ge k		;hipass filter
 ;use xm to mask each of the trig. arrays
 a=a*xm
 b=b*xm
 c=c*xm
 d=d*xm
 sc2db,xout,a,b,c,d
 return,xout
 endfunc
 ;=========================================================================
subr hipass,xin,xout,k
 ;makes a high pass filtered version of xin and puts it in xout
 ;the cutoff freq. # is k
 sc2d,xin,a,b,c,d
 nx=dimen(a,0)
 ny=dimen(a,1)
 maskprep,xm,nx,ny	;get array of freq. #'s
 xm=xm ge k		;hipass filter
 ;use xm to mask each of the trig. arrays
 a=a*xm
 b=b*xm
 c=c*xm
 d=d*xm
 sc2db,xout,a,b,c,d
 endsubr
 ;=========================================================================
subr hipassm,xx,k
 ;makes a high pass filtered version of xin and puts it in xout
 ;the cutoff freq. # is k
 ;this version doesn't save input, designed to use less memory
 sc2d,xx,a,b,c,d
 xx=0
 nx=dimen(a,0)
 ny=dimen(a,1)
 maskprep,xm,nx,ny	;get array of freq. #'s
 xm=xm ge k		;hipass filter
 ;use xm to mask each of the trig. arrays
 a=a*xm
 b=b*xm
 c=c*xm
 d=d*xm
 xm=0
 sc2db,xx,a,b,c,d
 endsubr
 ;=========================================================================
func lopass(xin,k)
 ;returns a low pass filtered version of xin
 ;the cutoff freq. # is k
 sc2d,xin,a,b,c,d
 nx=dimen(a,0)
 ny=dimen(a,1)
 maskprep,xm,nx,ny	;get array of freq. #'s
 xm=xm le k		;lopass filter
 ;use xm to mask each of the trig. arrays
 a=a*xm
 b=b*xm
 c=c*xm
 d=d*xm
 sc2db,xout,a,b,c,d
 return,xout
 endfunc
 ;=========================================================================
 subr lopass,xin,xout,k
 ;makes a low pass filtered version of xin and puts it in xout
 ;the cutoff freq. # is k
 sc2d,xin,a,b,c,d
 nx=dimen(a,0)
 ny=dimen(a,1)
 maskprep,xm,nx,ny	;get array of freq. #'s
 xm=xm le k		;lopass filter
 ;use xm to mask each of the trig. arrays
 a=a*xm
 b=b*xm
 c=c*xm
 d=d*xm
 sc2db,xout,a,b,c,d
 endsubr
 ;=========================================================================
subr lopassm,xx,k
 ;makes a lo pass filtered version of xin and puts it in xout
 ;the cutoff freq. # is k
 ;this version doesn't save input, designed to use less memory
 sc2d,xx,a,b,c,d
 xx=0
 nx=dimen(a,0)
 ny=dimen(a,1)
 maskprep,xm,nx,ny	;get array of freq. #'s
 xm=xm le k		;lopass filter
 ;use xm to mask each of the trig. arrays
 a=a*xm
 b=b*xm
 c=c*xm
 d=d*xm
 xm=0
 sc2db,xx,a,b,c,d
 endsubr
 ;=========================================================================
subr fft2d,x,cr,ci
 ;returns real and imag. parts of fourier transform coefficients
 ;4/9/90 - some errors in the imaginary part discovered
 !fftdp=0	;single precision seems adequate
 sc2d,x,a,b,d,c	;note reversal of c & d
 ;do the edges right
 ;since various terms were zero, this amounts to multiplying along the edges
 a(0,0)=2.*a(*,0)
 n=dimen(a,0)-1
 m=dimen(a,1)-1
 a(0,m)=2.*a(*,m)
 for i=0,m do { a(0,i)=2.*a(0,i)   a(n,i)=2.*a(n,i)  }
 ;note that the corners are multiplied by 4
 ;now also for c and d, b is 0 along all edges, c and d along some
 c(0,0)=2.*c(*,0)
 c(0,m)=2.*c(*,m)
 for i=0,m do { c(0,i)=2.*c(0,i)   c(n,i)=2.*c(n,i)  }
 d(0,0)=2.*d(*,0)
 d(0,m)=2.*d(*,m)
 for i=0,m do { d(0,i)=2.*d(0,i)   d(n,i)=2.*d(n,i)  }
 xq=a+b
 a=a-b
 b=-c-d
 c=d-c
 d=0
 cr=reverse(xq,1)
 cr=concat(cr,a(*,1:*))
 cr=cr(>1,>0)
 cr=concat(reverse(cr),cr(*,1:*))
 a=0
 xq=0
 ci=reverse(c,1)
 c=0
 ci=concat(ci,b(*,1:*))
 b=0
 ci=ci(>1,>0)
 ci=concat(-reverse(ci),ci(*,1:*))
 endsubr
 ;=========================================================================
subr filter,xout,xm,a,b,c,d
 ;generates xout after applying the filter xm to a,b,c,d
 sc2db,xout,a*xm,b*xm,c*xm,d*xm
 endsubr
 ;=========================================================================
subr unsharp,xin,xout,k
 ;performs unsharp masking with width of k
 xout=smooth(xin,k)
 xout=xout(>1,>0)
 xout=smooth(xout,k)
 xout=xin-xout(>1,>0)
 endsubr
 ;=========================================================================
func unsharp(xin,k)
 ;performs unsharp masking with width of k
 xout=smooth(xin,k)
 xout=xout(>1,>0)
 xout=smooth(xout,k)
 xout=xin-xout(>1,>0)
 return,xout
 endfunc
 ;=========================================================================
subr compower2,x,p
 ;generates a symmetrized power spectrum from the 2-D array x
 ;the result is a 1-D vector p
 sc2d,x,a,b,c,d
 a=a*a+b*b+c*c+d*d
 ;save memory
 b=0  c=0  d=0
 nx=dimen(a,0)  ny=dimen(a,1)
 maskprep,msk,nx,ny
 ;msk contains the freq. number of each element in the 2-D array a
 msk=msk+0.5	;add 0.5 so subsequent truncation is actually a rounding
 np=max(msk)+1
 p=fltarr(np)
 zero,p	;w will be used as the weight for each freq. (the number of times
 w=p	;it occurs in the 2-D array
 for i=0,nx-1 do for j=0,ny-1  do {
   ;load p and calculate w
   k=msk(i,j)
   w(k)+=1.0
   p(k)+=a(i,j)
 }
 ;now normalize p
 w=w>1.0	;if a w was zero, then that p will be also
 p=p/w
 endsubr
 ;=========================================================================
func smooth2,x,ns
 ;a slow way to boxcar smooth
 xs=smooth(x,ns)
 xs=smooth(xs(>1,>0),ns)
 xs=xs(>1,>0)
 return,xs
 endfunc
 ;=========================================================================
func radialsmooth(xin, rad)
 ;this assumes center at nx/2, ny/2 which is consistent with 2-D power arrays
 ;produced by power2dw
 ;make a version of xin that is radially smoothed, useful for subtracting, scaling,
 ;or thresholding
 if isarray(xin) eq 0 then { ty,'radialsmooth accepts only 2-D arrays' return, 0 }
 if num_dim(xin) ne 2 then { ty,'radialsmooth accepts only 2-D arrays' return, 0 }
 nx = dimen(xin,0)   ny = dimen(xin,1)
 rq = sqrt( (indgen(xin,0) - nx/2)^2 + (indgen(xin,1) - ny/2)^2)
 rq = rfix(rq)
 key_locations, rq, values, counts, indexsymarr
 ty,'total(counts) vs nx*ny ', fix(total(counts)), nx*ny
 np = num_elem(values)
 rad = fltarr(np)
 xout = fltarr(nx,ny)
 xout(nx/2,ny/2) = xin(nx/2,ny/2)
 for k = 0, np-1 do {
   inq = indexsymarr(k)
   xq = mean(xin(inq))
   xout(inq) = xq
   rad(k) = xq
 }
 return, xout
 endfunc
 ;=========================================================================
func radialmediansmooth(xin, rad)
 ;this assumes center at nx/2, ny/2 which is consistent with 2-D power arrays
 ;produced by power2dw
 ;this differs from radialsmooth, it does a median rather an average for each radii
 if isarray(xin) eq 0 then { ty,'radialmediansmooth accepts only 2-D arrays' return, 0 }
 if num_dim(xin) ne 2 then { ty,'radialmediansmooth accepts only 2-D arrays' return, 0 }
 nx = dimen(xin,0)   ny = dimen(xin,1)
 rq = sqrt( (indgen(xin,0) - nx/2)^2 + (indgen(xin,1) - ny/2)^2)
 rq = rfix(rq)
 key_locations, rq, values, counts, indexsymarr
 ty,'total(counts) vs nx*ny ', fix(total(counts)), nx*ny
 np = num_elem(values)
 rad = fltarr(np)
 xout = fltarr(nx,ny)
 xout(nx/2,ny/2) = xin(nx/2,ny/2)
 for k = 0, np-1 do {
   inq = indexsymarr(k)
   yq = sort(xin(inq))
   mid = num_elem(yq)/2
   xq = yq(mid)
   xout(inq) = xq
   rad(k) = xq
 }
 return, xout
 endfunc
 ;=========================================================================
func radialmedian4(xin, rad)
 ;this assumes center at nx/2, ny/2 which is consistent with 2-D power arrays
 ;produced by power2dw
 ;this differs from radialmediansmooth, it also uses medians but
 ;breaks up the circum. into 45 degree zones for better performance against rays
 if isarray(xin) eq 0 then { ty,'radialmedian4 accepts only 2-D arrays' return, 0 }
 if num_dim(xin) ne 2 then { ty,'radialmedian4 accepts only 2-D arrays' return, 0 }
 nx = dimen(xin,0)   ny = dimen(xin,1)
 rq = sqrt( (indgen(xin,0) - nx/2)^2 + (indgen(xin,1) - ny/2)^2)
 rq = rfix(rq)
 key_locations, rq, values, counts, indexsymarr
 ty,'total(counts) vs nx*ny ', fix(total(counts)), nx*ny
 np = num_elem(values)
 rad = fltarr(np)
 xout = fltarr(nx,ny)
 xout(nx/2,ny/2) = xin(nx/2,ny/2)
 for k = 0, np-1 do {
   inq = indexsymarr(k)
   if k le 5 then {
     yq = xin(inq)
     yq = sort(yq)
     mid = num_elem(yq)/2
     xq = yq(mid)
     xout(inq) = xq
     rad(k) = xq
   } else {
     x = inq%nx - nx/2
     y = inq/nx - ny/2
     th=atan2(x,y)
     ith = index(th)
     inq = inq(ith)
     yq = xin(inq)
     nq = num_elem(yq)
     ;break into 8 parts
     mid = nq/16
     rmean = 0
     for i=0,7 do {
       i0 = float(i*nq)*0.125 > 0
       i1 = float((i+1)*nq - 1)*0.125 < (nq-1)
       zq = sort(yq(i0:i1))
       xq = zq(mid)
       xout(inq(i0:i1)) = xq
       rmean += xq
     }
     rad(k) = 0.125 * rmean
   }
 }
 return, xout
 endfunc
 ;=========================================================================
func radialthetasmooth(xin, sm)
 ;this assumes center at nx/2, ny/2 which is consistent with 2-D power arrays
 ;produced by power2dw
 ;this differs from radialmediansmooth, it also uses medians but
 ;breaks up the circum. into 45 degree zones for better performance aginst rays
 if isarray(xin) eq 0 then { ty,'radialthetasmooth accepts only 2-D arrays' return, 0 }
 if num_dim(xin) ne 2 then { ty,'radialthetasmooth accepts only 2-D arrays' return, 0 }
 nx = dimen(xin,0)   ny = dimen(xin,1)
 rq = sqrt( (indgen(xin,0) - nx/2)^2 + (indgen(xin,1) - ny/2)^2)
 rq = rfix(rq)
 key_locations, rq, values, counts, indexsymarr
 ty,'total(counts) vs nx*ny ', fix(total(counts)), nx*ny
 np = num_elem(values)
 xout = fltarr(nx,ny)
 xout(nx/2,ny/2) = xin(nx/2,ny/2)
 for k = 0, np-1 do {
   inq = indexsymarr(k)
   if k le 5 then {
     yq = xin(inq)
     yq = sort(yq)
     mid = num_elem(yq)/2
     xq = yq(mid)
     xout(inq) = xq
   } else {
     x = inq%nx - nx/2
     y = inq/nx - ny/2
     th=atan2(x,y)
     ith = index(th)
     inq = inq(ith)
     yq = xin(inq)
     yq = gsmooth(yq, sm)
     xout(inq) = yq
   }
 }
 return, xout
 endfunc
 ;=========================================================================