;WRITTEN: James Mason 2008-2009 (jmason86@gmail.com)
;MODIFIED: James Mason 1/21/2010
;
;
; Given a single AR id #:
; - finds all MDI images containing that AR
; - copies from MDI database and renames fits files
; - uses pe and fastrack modules to make projection and extract region

@/auto/home2/jmason86/Forecast/pros/sel_AR.pro
@/auto/home2/jmason86/Forecast/pros/mdiquality.pro
@/auto/home2/jmason86/Forecast/pros/bzcorrect.pro
@/auto/home2/jmason86/Forecast/pros/mapout.pro
@/auto/home2/jmason86/Forecast/pros/ARinfo.pro

pro findar, no_copy=no_copy, out_typ=out_typ, dat_out=dat_out, gif=gif, ps=ps, loud=loud, typ=typ, lat=lat, lon=lon, area=area, lng_ext=lng_ext, lbl=lbl, src=src, z=z, long_cut=long_cut, fix_dir=fix_dir, yn=yn, nar=nar, sorc=sorc, seldate=seldate, debug=debug, sav_fits=sav_fits

IF n_elements(debug) eq 0 then debug = 0
; create and set location of working directory
IF debug then goto, debuger
working='/auto/tmp21/jmason86/findAR'
cd, working
;spawn,'pwd > dir'
;readcol, 'dir', working, format='a',/silent
;working = working(0)
;spawn, 'rm dir'


;determine NAR to run on
if n_elements(noaa) eq 0 then begin
print, "Select ARs by: 1) NOAA AR#'s, 2) Date Range (>1 week ago), 3) Most Recent ARs 4) Range of ARs"
choice=0 & read, choice
if choice eq 1 then goto, by_AR & if choice eq 2 then goto, by_date & if choice eq 3 then goto, by_yester & if choice eq 4 then goto, by_rng

by_AR:
iii=0
print, 'How many ARs?:'
totar = 0 & read, totar
noaa_all = intarr(totar)
print, 'Input NOAA AR#:'
read, noaa_all
data_loc = STRING('/mag')
goto, done

by_date:
iii=0
print, 'Input start date. Format = yyyy-MON-dd'
seldate = STRING('') & read, seldate
print, 'Input end date. Format = yyyy-MON-dd'
seldate2 = STRING('') & read, seldate2
data_loc = STRING('/mag')
narset = get_nar(seldate,seldate2)
noaa_all = narset.noaa

goto, done

by_yester:
print, 'Recent data compiling...'
seldate = STRING('')
data_loc = STRING('/synoptic/eof')
narset = get_nar(seldate)
noaa_all = narset.noaa
goto, done

by_rng:
iii=0
print, 'Input first NOAA AR#:'
firar = 0 & read, firar
print, 'Input last NOAA AR#;'
lasar = 0 & read, lasar
totar = lasar - firar + 1
noaa_all = intarr(totar)
FOR alpha = 0, totar-1 DO BEGIN
  noaa_all(alpha) = firar + alpha
ENDFOR
data_loc = '/mag'
goto, done 


debuger: 
working = '/auto/home2/jmason86/solignis/solignis/data'
yn = 'n'
noaa_all=7999
data_loc = '/mag'
cd, working

done:
endif

;NAR loop
priorNAR = intarr(n_elements(noaa_all)) & priorNAR(*) = !VALUES.F_NAN & ii=0
print, 'NOAA Active regions to run:', noaa_all
last = (n_elements(noaa_all)-1)
REPNAR:
WHILE iii LE last do begin
prior_nar_check = where(priorNAR EQ noaa_all(iii))
IF prior_nar_check EQ [-1] THEN noaa = noaa_all(iii) ELSE BEGIN 
  iii=iii+1
  goto, REPNAR
ENDELSE
IF noaa_all(iii) EQ [0] OR noaa_all(iii) EQ [8253] OR noaa_all(iii) EQ [8255] THEN BEGIN
  iii=iii+1
  goto, REPNAR
ENDIF

priorNAR(ii) = noaa_all(iii)


print,'Program running on AR:', noaa


defaults:
; set call parameter defaults
if n_elements(no_copy) eq 0 then no_copy=0
if n_elements(out_typ) eq 0 then out_typ=1
if n_elements(dat_out) eq 0 then dat_out=1
if n_elements(gif) eq 0 then gif=1
if n_elements(ps) eq 0 then ps=0
if n_elements(loud) eq 0 then loud=0
if n_elements(typ) eq 0 then typ=0
if n_elements(z) eq 0 then z=0
if n_elements(long_cut) eq 0 then long_cut=40
if typ eq 1 and n_elements(src) eq 0 then no_copy=1
if gif eq 0 and ps eq 0 then begin & gif=1 & endif

lbl=strcompress(string(fix(noaa)),/remove_all)
spawn, 'rm -rf '+lbl
spawn, 'mkdir '+lbl


; create directories to store temporary files in
if (no_copy eq 0) then spawn, 'mkdir '+lbl+'/fits'
spawn, 'mkdir out'
if n_elements(src) ne 0 then begin & spawn, 'cp '+src+'/*.fits ./'+lbl+'/fits/' & no_copy=1 & endif

; find the range of mdi images that the AR appears in
; loud turns on feedback messages in sel_AR program
; long_cut sets the boundaries of longitude the AR must be inside
if (typ eq 0) then begin 
sel_AR, noaa, srt_dt, subsq, edpt, errdat, dims, loud=loud, long_cut=long_cut
  if errdat eq 1 then begin
    error = 'Not contained in MDI data'
    goto, fin2
  endif
endif else if (typ eq 1) then begin
  dims=fltarr(6)
  if n_elements(lat) eq 0 then begin & print, "latitude:" & lat=0.0 & read, lat & endif
  dims(2)=lat
  if n_elements(lon) eq 0 then begin & print, "carrington longitude:" & lon=0.0 & read, lon & endif
  dims(3)=lon
  if n_elements(area) eq 0 then begin & print, "area (in 1e-6 * FD):" & area=0.0 & read, area & endif
  dims(1)=area
  if n_elements(lng_ext) eq 0 then begin & print, "longitudinal extent in degrees:" & lng_ext=0.0 & read, lng_ext & endif
  dims(0)=lng_ext 
  if n_elements(lbl) eq 0 then begin & print, "output folder label" & lbl='' & read, lbl & endif
  errdat=0
endif
if errdat eq 1 then goto, fin2

; using time_index module, converts the start date found to a SOHO era day
if (no_copy eq 0) then begin
dt=strcompress(string(srt_dt(0)),/remove_all)+'.'+strcompress(string(srt_dt(1)),/remove_all)+'.'+strcompress(string(srt_dt(2)),/remove_all)
spawn, 'rm -f tmp'
spawn, 'time_index -d time='+dt+' > tmp'
close, 5
openr,5,'tmp'
srt_num=0 & readf,5,srt_num
close,5
num=srt_num
idx=0 & ctr=1 & tst=0
subsq=[subsq, 0, 0]
if subsq(0) eq 0 then begin & idx=2 & endif

; copies mdi fits files in whole day blocks from database to temporary folder src
repeat begin
  number=strcompress(string(num),/remove_all)
  spawn, 'ls /synoptic/eof/fd_M_96m_01d.00'+number+'/fd_M_96m_01d.'+number+'.00*.fits >'+lbl+'dlist_mvp'
  log_name=+lbl+'dlist_mvp'
  close, 50 & openr,50,log_name
	  while not eof(50) do begin
	    data_loc=STRING('/synoptic/eof')
          endwhile
  close, 50
  print, 'cp '+data_loc+'/fd_M_96m_01d.00'+number+'/fd_M_96m_01d.'+number+'.00*.fits ./'+lbl+'/fits/'
  spawn, 'cp '+data_loc+'/fd_M_96m_01d.00'+number+'/fd_M_96m_01d.'+number+'.00*.fits ./'+lbl+'/fits/'
  if ctr lt subsq(idx) then begin & ctr=ctr+1 & num=num+1 & endif else begin 
  ctr=1 & idx=idx+1 
  if (n_elements(subsq) gt idx) then num=num+subsq(idx) & idx=idx+1 & endelse
  if n_elements(subsq) le idx then tst=1 else if subsq(idx) eq 0 then tst=1
STOP
endrep until tst eq 1

; cuts out files where AR is outside boundaries in first and last days' files
id=-1 & if edpt(0) gt 0 then repeat begin
id=id+1
if id lt 10 then idn='0'+strcompress(string(id),/remove_all) else idn=strcompress(string(id),/remove_all)
;	print, 'rm -f ./'+lbl+'/fits/fd_M_96m_01d.'+strcompress(string(srt_num),/remove_all)+'.00'+idn+'.fits'
;spawn, 'rm -f ./'+lbl+'/fits/fd_M_96m_01d.'+strcompress(string(srt_num),/remove_all)+'.00'+idn+'.fits'
endrep until id eq edpt(0)-1

id=15 & if edpt(1) lt 14 then repeat begin
id=id-1
if id lt 10 then idn='0'+strcompress(string(id),/remove_all) else idn=strcompress(string(id),/remove_all)
print, 'rm -f ./'+lbl+'/fits/fd_M_96m_01d.'+strcompress(string(num),/remove_all)+'.00'+idn+'.fits'
spawn, 'rm -f ./'+lbl+'/fits/fd_M_96m_01d.'+strcompress(string(num),/remove_all)+'.00'+idn+'.fits'
endrep until id eq edpt(1)+1 
endif

;if no_src eq 1 then begin & del_stuff=1 & goto, no_files & endif

; makes directory to store output in
spawn, 'mkdir '+lbl

;Determine good data/drop bad data
        ;list fits files to get mdi_id from file name
spawn, 'ls ./'+lbl+'/fits/*.fits >> ./'+lbl+'/fits/list'
openr, 6, './'+lbl+'/fits/list'
array = strarr(100000) & ar=0L & ar_end=ar
while not eof(6) do begin
file  = '' & readf, 6, file
array(ar) = file

fits = readfits(array(ar),hd, /silent)
DATAFILE = strarr(27)
d=fxpar(hd,'DATAFILE')
mdi_id = strmid(d, 13, 4)
fits_id = strmid(d, 18, 4)
t=fxpar(hd,'T_REC')
time = strmid(t,0,23)

        ;determine mdi quality
record='/auto/PAS/D136007/mdi_log/lev1.5/fd_M_96m_01d/00'+mdi_id+'.record.rdb'
info=strarr(17)
mag_qual=strarr(15) & mag_vld=strarr(15)
mag_qual(0)='NoFile'

openr, 1, record, error=errout
if errout eq 0 then begin
readf, 1, info
quality = 0L & qualvld = 0L & mag_qual=0L & mag_vld=0L
pos1=strpos(info(0),'QUALITY')
val=strpos(info(0),'QUAL_VLD')
timepos=strpos(info(0),time)
FOR j=2, 16 do begin
time2=strmid(info(j), timepos, 23)
if time EQ time2 then begin
    mag_qual=strmid(info(j), pos1, 10)
    reads, mag_qual, quality, format='(%"%x")'
    mag_vld=strmid(info(j), val, 10)
    reads, mag_vld, qualvld, format='(%"%x")'
endif
ENDFOR ; j loop
maskstring = '0x402c0c04' & qualmask = 0L
reads, maskstring, qualmask, format='(%"%x")'
   endif

bad = 0
if ((quality AND qualmask AND qualvld) NE 0) then bad = 1 
	if bad eq 1 then begin
	  spawn, 'mkdir ./'+lbl+'/fits/bad_data' 
	  spawn, 'mv ./'+lbl+'/fits/fd_M_96m_01d.'+mdi_id+'.'+fits_id+'.fits ./'+lbl+'/fits/bad_data'
	  print, 'Removing bad data: ', d
	endif
close,1
ar=ar+1L & ar_end=ar-1L
endwhile & close, 6
spawn, 'rm ./'+lbl+'/fits/list'

;end determine good data/drop bad data

; gets list of files in src folder
no_src=1
spawn,'\rm '+lbl+'dlist_mvp'
spawn,'ls ./'+lbl+'/fits/*.* > '+lbl+'dlist_mvp'
close, 5
openr,5,+lbl+'dlist_mvp'
dfname = strarr(10000) & dfp=0L & dfp_end=dfp
while not eof(5) do begin
no_src=0
data_file='' & readf, 5, data_file
dfname(dfp) = data_file
dfp=dfp+1L & dfp_end=dfp-1L
endwhile & close,5
                                                                                                                                                             
del_stuff = 0 & if no_src eq 1 then del_stuff=1 

; starts output of analysis data

; analyses the size of the AR over its existence and picks a good size
; for the region of analysis to extract, maintains consistent sample size
; over analysis run

data = readfits(dfname(0), hd, /silent)
xscale=fxpar(hd,'XSCALE')
yscale=fxpar(hd,'YSCALE')

if (n_elements(where(dims))/5-1) gt 0 then if n_elements(dims(0,*)) gt 1 then dims=dims(*,0:n_elements(where(dims))/5-1)

sub=0 & if n_elements(dims(0,*)) gt 1 then long_ext=max(dims(0,0:n_elements(dims(0,*))-2), sub) else lo
ng_ext=dims(0)
area=dims(1,sub) & lat=dims(2,sub)
width=(long_ext>1) * 8.55 * cos(lat*3.14159/180) * xscale/2 + 50

if area le 40 then begin
  err_dat = 1
  print, "skipping: size too small"
  error = '' & error = 'size too small'
endif else err_dat = 0
if err_dat eq 1 then  goto, fin2

area=(area > 10)*xscale*yscale
height=area/width+100
day_ctr=intarr(2) & day_ctr(0)=-1 & day_ctr(1)=0

; begins loops that processes each fits file separately
FOR ip=0L, dfp_end do begin
; reads fits file
data = readfits(dfname(ip), hd, /silent)
; reads date and time from header
d=fxpar(hd,'T_REC')
if day_ctr(1) ne fix(strmid(d,8,2)) then begin & day_ctr(0)=day_ctr(0)+1 & day_ctr(1)=fix(strmid(d,8,2)) & endif
dt=strmid(d, 0, 4)+'.'+strmid(d,5,2)+'.'+strmid(d,8,2)
tm=strmid(d,11,12)
xscale=fxpar(hd,'XSCALE')
yscale=fxpar(hd,'YSCALE')
tmpname=fxpar(hd,'DATAFILE')


lat=median(dims(2,0:max(where(dims(2,*)))))
lon=median(dims(3,0:max(where(dims(3,*)))))
; assigns file number
if ip lt 10 then nm='000'+strcompress(string(ip),/remove_all) else if ip lt 100 then nm='00'+strcompress(string(ip),/remove_all) else if ip lt 1000 then nm='0'+strcompress(string(ip),/remove_all) else nm=strcompress(string(ip),/remove_all)

; corrects Bz values of pixels in image using helioscopic angle
bzcorrect, +lbl+'/fits/'+tmpname

spawn, 'mkdir '+lbl+'/sub/'

; creates tst.map file for running pe module
mapout, fn=+lbl+'tst', in= +working+'/'+lbl+'/fits/'+tmpname, out= +working+'/'+lbl+'/sub/', rows=fix(height), cols=fix(width), lat=lat, lon=lon


; runs project module using pe and tst.map file
spawn, 'pe map='+lbl+'tst.map'

; renames pe output so reiteration will not overwrite
cd, working
spawn, 'mv ./'+lbl+'/sub/0000.fits ./'+lbl+'/sub/NOAA'+lbl+'_'+nm+'_'+dt+'_'+tm+'.fits'


IF err_dat EQ 0 THEN ARinfo, fn=lbl, rows=fix(height), cols=fix(width), lat=lat, lon=lon

ENDFOR ;ip loop

fin2:
if n_elements(err_dat) eq 0 then err_dat = 0 & if err_dat eq 1 then begin 
  cd, working
  spawn, 'rm -rf '+lbl 
  close, 10 & if iii eq 0 then openw, 10, 'skipped_regions' else $
                               openw, 10, 'skipped_regions', /append
  printf, 10, lbl, ' ', error & print, lbl, ' ', error
  close, 10
endif

;cleanup
cd, working
spawn, 'mkdir junk' & spawn, 'mv *.param *.gif *tmp *mvp junk'
spawn, 'rm -rf junk'
spawn, 'rm '+working+'/'+lbl+'tst.map' & spawn, 'rm -rf out'
IF not keyword_set(sav_fits) THEN spawn, 'rm -rf '+lbl+'/fits'
spawn,'mv '+lbl+' NOAA'+lbl

iii=iii+1 & ii=ii+1
ENDWHILE   ;iii loop 


if errdat eq 1 then print, 'No AR data found for '+lbl+' in NOAA database' else print, '********************Program normal completion********************'


END