;this assumes that you have run HMIfilters.pro, and somehow saved the output.
PRO APPLY_HMI_PROFILES,stokesfile,filterfile,nx,ny,nw,dlam,outfile,yoff=yoff

; stokesfile = file of stokes spectra, each spatial point having [nw,4] size spectra
; filterfile = file where hmi filter information is saved for data having same nx,ny as
;	       stokesfile
; outfile = where to write the nx,ny,nfilt,4 data
; nx,ny size of model data and filter results.
; nw is number of wavelengths in spectra & filter data
; dw is spectral sampling interval
; nfilt for hmi hardwired to 6!!! 
; yoff is starting y (row) for the split-up 4096^2 files.  use ny=real ny (e.g., 512).

nfilt=6
yo=0L
if (keyword_set(yoff)) then yo=long(yoff)

openr,lun1,stokesfile,/get_lun
openr,lun2,filterfile,/get_lun
openw,lun3,outfile,/get_lun

a=assoc(lun1,fltarr(nw,4)) ; 4 = stokes....
b=assoc(lun2,fltarr(nw,nfilt))
hmidat=fltarr(nfilt,4)

progress,0.,/reset,label='HMI filters: ',small=5.,big=10.
nprogress = 0.0d0
fprogress = double(nx)*double(ny)

for j=0L,ny-1 do begin ; work up the rows
   for i=0L,nx-1 do begin ; work up the columns, x is fast-axis
      inspect=a(j*nx + i)
      ;inspect=a(i*nx + j)
      ;filt=b(j*nx + i)
      filt=b(i*nx + j+yo)
      for k=0,nfilt-1 do begin  ; loop through filters
         hmidat(k,0)=TOTAL(filt[*,k]*inspect[*,0])*dlam
         hmidat(k,1)=TOTAL(filt[*,k]*inspect[*,1])*dlam
         hmidat(k,2)=TOTAL(filt[*,k]*inspect[*,2])*dlam
         hmidat(k,3)=TOTAL(filt[*,k]*inspect[*,3])*dlam
      endfor ; filter loop
      writeu,lun3,hmidat
      nprogress = temporary(nprogress) + 1.0d0
      progress,100.0*(nprogress/fprogress)
   endfor ; x-dir loop
endfor ; y-dir loop

close,lun1 & free_lun,lun1
close,lun2 & free_lun,lun2
close,lun3 & free_lun,lun3

progress,100.0,/last


end