; Find HMI magnetic field image at closest to XRT
; IDL> .r find_hmi_mag_qsbp
;
;--- read XRT qsbp info structure ---
in_dir='/disk/yla/takeda/REU_2026/qsbp_info/'
fname='qsbp_info_20230219'
restgen,file=in_dir+fname,qsbp_info
yyyymmdd=strmid(fname,10,8)
fov_xyr=qsbp_info.fov_xyr
bp_xyr=qsbp_info.bp_xyr
;
;--- read XRT images ---
file_ap=qsbp_info.file_ap
mreadfits,file_ap,index,data
;
;*** shift Al_poly image: Sun-center to image center ***
sftx1=round(index.xcen/index.xscale)
sfty1=round(index.ycen/index.yscale)
dd2_ap=shift(data,sftx1,sfty1)
;
;*** rotate Al_poly images: Sun's North is up ***
ang1=index.crota2*(-1)
magfact1=1.0 & cx=512 & cy=512
dd3_ap=rot(dd2_ap,ang1,magfact1,cx,cy,/pivot,/interp)
;
;stop
; goto,HIST
;--- find HMI image  ---
hmi_dir='/disk/data/SDO/SG/hmi.M_720s/'
dt0=3600.
r_val=1.0  ; use same value as XRT
;r_val=0.8
b_cutoff=10.

h_scl0=0.60
cc0=(h_scl0*725e5)^2.   ; Gauss --> Mx

; 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
 y_str=strmid(index(i).date_obs,0,4)
 m_str=strmid(index(i).date_obs,5,2)
 d_str=strmid(index(i).date_obs,8,2)

 h_dir=hmi_dir+y_str+'/'+m_str+'/'+d_str+'/'
 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))
 h_da(*,*,i)=h_data

endfor
;
;--- shrink HMI to XRT size ---
ang=0.0 & hcx=2048 & hcy=2048
magfact=h_index.cdelt1/index.cdelt1
h_da2=rot(h_data,ang,magfact,hcx,hcy,/pivot,/interp)
h_da3=h_da2(hcx-512:hcx+512-1,hcy-512:hcy+512-1)

wdef,1,1024,1024
;
set_plot,'z'
device,set_resolution=[1024,1024]
loadct,0
tvlct,rr,gg,bb,/get
;
erase
tvscl,h_da3>(-50)<50
;
xyouts,10,10,/dev,'SDO/HMI mag '+anytim(h_in.date_obs,/ecs)
snap0=tvrd()
;
erase
tvscl,h_da3>(-50)<50
sftx=fov_xyr(0)-512
sfty=fov_xyr(1)-512
nn_bp=n_elements(bp_xyr(0,*))
;
draw_circle,fov_xyr(0),fov_xyr(1),fov_xyr(2),/dev
for i=0,nn_bp-1 do $
  draw_circle,bp_xyr(0,i)-sftx,bp_xyr(1,i)-sfty,bp_xyr(2,i),/dev
;
xyouts,10,10,/dev,'SDO/HMI mag '+anytim(h_in.date_obs,/ecs)
snap1=tvrd()
;
erase
tvscl,safe_log10(dd3_ap)
xyouts,10,10,/dev,'Hinode/XRT Op/Apol '+anytim(index.date_obs,/ecs)
snap2=tvrd()
;
draw_circle,fov_xyr(0),fov_xyr(1),fov_xyr(2),/dev
for i=0,nn_bp-1 do $
  draw_circle,bp_xyr(0,i),bp_xyr(1,i),bp_xyr(2,i),/dev
;
snap3=tvrd()
;
device,/close & set_plot,'x'
tvscl,snap1
;
wdef,0,1024,1024
tvscl,snap3
;
;*** for comparison ***
;  stepper,[[[snap3]],[[snap1]]]
comparr=bytarr(1024,1024,2)
comparr(*,*,0)=snap1
comparr(*,*,1)=snap3
stepper,comparr
;
;stop,'*** write PNG?'
;
;write_png,'./hmi0_'+yyyymmdd+'.png',snap0,rr,gg,bb
;write_png,'./hmi_'+yyyymmdd+'.png',snap1,rr,gg,bb
;write_png,'./xrt0_'+yyyymmdd+'.png',snap2,rr,gg,bb
;write_png,'./xrt_'+yyyymmdd+'.png',snap3,rr,gg,bb
;
END
