; Find HMI magnetic field image at closest to XRT and
;  - correct projection effect
;  - convert flux density (gauss) to flux (Mx) with the correction of
;     seasonal variation of pixel size (in cm^2) due to the change of 
;     Sun-Eatrh distance.
;  - mask out pixels with B cut-off less than 10 Mx (noise level).
;  - then get total unsigned magnetic flux within FOV, QS and BPs.
;  - save HMI file name and measured flux.
;
; IDL> .r calc_hmi_mag_unsig
;     (to read Nathaniel's save file, need to read XRT files via
;       /www/HINODE/XRT/SCIA. use Lewis, Filament, Illyrica, etc.)
; 
; 2026/06/20 :  initial version. 
; 2026/07/23 :  added Npix to the output structure.
;
;--- read XRT qsbp info structure ---
in_dir='/disk/yla/takeda/REU_2026/qsbp_info/'
; fname='qsbp_info_20230219'  ; sample case for non-Min.
; fname='qsbp_info_20191219'  ; sample case for SolMin.
;-------------------------------
fname='qsbp_info_20111227'
fname='qsbp_info_20130201'
;     fname='qsbp_info_20150723'
;fname='qsbp_info_20171229'
;-------------------------------
;fname='qsbp_info_20180709'
;     fname='qsbp_info_20181201'
;fname='qsbp_info_20190519'
;fname='qsbp_info_20190917'
;fname='qsbp_info_20191219'
;fname='qsbp_info_20200212'
;fname='qsbp_info_20200513'
;fname='qsbp_info_20200911'
;fname='qsbp_info_20210108'
;fname='qsbp_info_20210207'
;-------------------------------
;fname='qsbp_info_20220624'
;fname='qsbp_info_20230219'
;fname='qsbp_info_20241116'
;-------------------------------
;
restgen,file=in_dir+fname,qsbp_info
yyyymmdd=strmid(fname,10,8)
file_ap=qsbp_info.file_ap
fov_xyr=qsbp_info.fov_xyr
bp_xyr=qsbp_info.bp_xyr
nn_bps=n_elements(bp_xyr(0,*))
;
;--- read Al_poly image to get header info, xcen ycen, etc.
mreadfits,file_ap,index,/nodata
;
;--- convert from (XRT image pix) to (arcsec from sun center) ---
; [Note] XRT image was roll angle corrected and shifted to image 
;     center when determining fov_xyr and bp_xyr.
;
fov_xyr_arcs=fov_xyr
bp_xyr_arcs=bp_xyr
;
xrt_cx=512      ; -index.xcen/index.xscale
xrt_cy=512      ; -index.ycen/index.yscale
fov_xyr_arcs(0)=(fov_xyr(0)-xrt_cx)*index.xscale
fov_xyr_arcs(1)=(fov_xyr(1)-xrt_cy)*index.yscale
fov_xyr_arcs(2)=fov_xyr(2)*index.xscale
;
for i=0,nn_bps-1 do begin
   bp_xyr_arcs(0,i)=(bp_xyr(0,i)-xrt_cx)*index.xscale
   bp_xyr_arcs(1,i)=(bp_xyr(1,i)-xrt_cy)*index.yscale
   bp_xyr_arcs(2,i)=bp_xyr(2,i)*index.xscale
endfor
;
;******** HMI image processing *********
hmi_dir='/disk/data/SDO/SG/hmi.M_720s/'
dt0=3600.
r_val=1.0
;r_val=0.8
b_cutoff=10.
;b_cutoff=40.

h_scl0=0.60   ; h_index.cdelt1   (arcsec/pix)

; Prepare for mask
sz0=4096
tmp_xx=findgen(sz0) # (fltarr(sz0)+1)
tmp_yy=(fltarr(sz0)+1) # findgen(sz0)
dd_arr=sqrt((tmp_xx-(sz0-1)/2.)^2.+(tmp_yy-(sz0-1)/2.)^2.)*h_scl0

nn0=n_elements(index)
tot_magflx=dblarr(nn0)
h_da=fltarr(sz0,sz0,nn0)
h_ff=strarr(nn0)

for i=0,nn0-1 do begin
 yyyy=strmid(yyyymmdd,0,4)
 mm=strmid(yyyymmdd,4,2)
 dd=strmid(yyyymmdd,6,2)
;
 h_dir=hmi_dir+yyyy+'/'+mm+'/'+dd+'/'
 h_files=find_files('*.fits',h_dir)
 read_sdo,h_files,h_in

 ss=where(h_in.quality eq 0 or h_in.quality eq 128 or h_in.quality eq 1024 or $
    h_in.quality eq 4096 or h_in.quality eq 65536 or h_in.quality eq 66560,nn)

  if nn eq 0 then begin
   print,'No HMI data Matched'
   stop
  endif

 h_in=h_in(ss)
 h_files=h_files(ss)

 x_time=anytim(index(i).date_obs)
 h_times=anytim(h_in.date_obs)

 tmp_min=min(abs(h_times-x_time),tmp_ss)
  if tmp_min gt dt0 then begin
   print,'No HMI data Matched'
   stop
  endif

 read_sdo,h_files(tmp_ss(0)),in0,da0,/use_share
 hmi_prep,in0,da0,h_index,h_data
 h_ff(i)=h_files(tmp_ss(0))

 ; For Correction (projection effect & annual distance variation)
 rr0=in0.rsun_obs
 cos_v=sqrt((rr0^2.-dd_arr^2.)>0)/rr0
 tmp_ss=where(cos_v gt 0.01)
 mod_arr=fltarr(sz0,sz0)
 mod_arr(tmp_ss)=(1./cos_v(tmp_ss))
 cc0=(h_index.cdelt1/3600.*!dtor*h_index.dsun_obs*100.)^2.

 nan_ss=where(finite(h_data,/nan),nan_nn)
 if nan_nn ne 0 then h_data(nan_ss)=0
 tmp_ss=where(dd_arr gt (rr0*r_val),tmp_nn)
 if tmp_nn ne 0 then h_data(tmp_ss)=0

 h_data=h_data*mod_arr
 bc_mask=abs(h_data) ge b_cutoff

 ;h_da(*,*,i)=h_data

endfor

;
; stop
;**************************************************************
fov_xyr_hpix=fov_xyr
bp_xyr_hpix=bp_xyr
;
fov_xyr_hpix(0)=2047.5+fov_xyr_arcs(0)/h_index.cdelt1
fov_xyr_hpix(1)=2047.5+fov_xyr_arcs(1)/h_index.cdelt2
fov_xyr_hpix(2)=fov_xyr_arcs(2)/h_index.cdelt2
;
for i=0,nn_bps-1 do begin
   bp_xyr_hpix(0,i)=2047.5+bp_xyr_arcs(0,i)/h_index.cdelt1
   bp_xyr_hpix(1,i)=2047.5+bp_xyr_arcs(1,i)/h_index.cdelt1
   bp_xyr_hpix(2,i)=bp_xyr_arcs(2,i)/h_index.cdelt1
endfor
;
wdef,0,1024,1024
tvscl,rebin(h_data,1024,1024)>(-50)<50
draw_circle,fov_xyr_hpix(0)/4,fov_xyr_hpix(1)/4,fov_xyr_hpix(2)/4,/dev
for i=0,nn_bps-1 do $
   draw_circle,bp_xyr_hpix(0,i)/4,bp_xyr_hpix(1,i)/4,bp_xyr_hpix(2,i)/4,/dev
;
help,fov_xyr_hpix,bp_xyr_hpix
;
;***** make a mask and get total unsigned magnetic field *****
mask_fov=blank_circle(4096,4096, $
          fov_xyr_hpix(0),fov_xyr_hpix(1),fov_xyr_hpix(2),/image)*1.0
usmag_fov=total(abs(h_data)*mod_arr*mask_fov*bc_mask)*cc0
ss=where(mask_fov*bc_mask eq 1,nss)
npix_fov=nss
;usmag_fov=total(abs(h_data)*mask_fov)
;
usmag_bps=dblarr(nn_bps)
npix_bps=lonarr(nn_bps)
mask_qs0=make_array(4096,4096,/byte,value=1)
for i=0,nn_bps-1 do begin
  mask_bp=blank_circle(4096,4096, $
              bp_xyr_hpix(0,i),bp_xyr_hpix(1,i),bp_xyr_hpix(2,i),/image) 
  usmag_bps(i)=total(abs(h_data)*mod_arr*mask_bp*bc_mask)*cc0
  ss=where(mask_bp*bc_mask eq 1,nss)
  npix_bps(i)=nss
  mask_qs0=mask_qs0*(1-mask_bp)
endfor
;
mask_qs=mask_qs0*mask_fov
usmag_qs=total(abs(h_data)*mod_arr*mask_qs*bc_mask)*cc0
ss=where(mask_qs*bc_mask eq 1,nss)
npix_qs=nss
; usmag_qs0=usmag_fov-total(usmag_bps)
;
file_hmi=h_ff
;
help,file_hmi
help,usmag_fov,usmag_qs,usmag_bps,npix_qs,npix_bps
;
;--- make final info structure ---
usmag_hmi=create_struct('file_hmi',file_hmi, $
              'fov_xyr_hpix',fov_xyr_hpix,'bp_xyr_hpix',bp_xyr_hpix, $
              'usmag_fov',usmag_fov, 'usmag_qs',usmag_qs, $
              'usmag_bps',usmag_bps, 'npix_fov',npix_fov, $
              'npix_qs',npix_qs, 'npix_bps',npix_bps)
help,usmag_hmi,/str
;
yn='n'
read,prompt='*** save info_structure? [y/n] ',yn
if yn eq 'y' then begin $
  outdir='/disk/yla/takeda/REU_2026/qsbp_usmag/'
  fname='usmag_hmi_'+yyyymmdd
  text='Created by CALC_HMI_MAG_UNSIG.PRO, '+ $
    'with the b_cutoff value of '+strtrim(string(b_cutoff),2)+'. ' + $
    'The structure contains the info of file_hmi(HMI filename), '+ $
    'fov_xyr_hpix(FOV center_x, center_y, radius in HMI pix), '+ $
    'bp_xyr_hpix(BP center_x, center_y, radius in HMI pix), '+ $
    'usmag_fov(Unsigned magnetic flux in unit of Max for FOV data), '+ $
    'usmag_qs(Unsigned magneric flux for QS data), '+ $
    'usmag_bps(Unsigned magnetic flux for BPs data), '+ $
    'npix_fov(Number of non-zero pix within FOV area), '+ $
    'npix_qs(Number of non-zero pix within QS area), '+ $
    'npix_bps(Number of non-zero pix within BPs areas).'
  name=['usmag_hmi']
  savegen,usmag_hmi,file=outdir+fname,text=text,name=name
endif
;
FIN:
end