PRO plotj, lbl, loc, dt
print, 'Plotting.'
spawn, 'mkdir plots'
spawn, 'rm plots/*'

;match MDI date and GOES date
year = strmid(dt,0,4)
close, 1 & openr, 1, +loc+'/GOES_data/xray'+year
openw, 2,'calc/Flares.dat'
readcol, +loc+'/GOES_data/xray'+year, long_year, st_time, end_time, max_time,lat, let_class, num_class, format = 'a,a,a,a,a,a,a', /silent
file_length = n_elements(long_year)-1
    ; loop to compare GOES date to MDI date
for i = 0, file_length do begin
if strmid(long_year(i),5,2) eq strmid(dt,2,2) and strmid(long_year(i),7,2) eq strmid(dt,5,2) and strmid(long_year(i),9,2) eq strmid(dt,8,2) then begin ; if there is a flare on the same date of the current magnetogram
    ; fix GOES data column inconsistency
  if strlen(lat(i)) le 3 then begin 
    num_class(i) = let_class(i)
    let_class(i) = lat(i) 
  endif
    ;GOES date 
  if strmid(long_year(i),5,1) eq '9' then $ 
    GOES_year = '19' + strmid(long_year(i), 5,2) else $
    GOES_year = '20' + strmid(long_year(i),5,2)
    GOES_month = strmid(long_year(i),7,2) & GOES_day = strmid(long_year(i),9,2)  
  GOES_date = GOES_year + '-' + GOES_month + '-' + GOES_day
  if let_class(i) ge 'M' then $ ; only include flares class M or higher   
  printf,2, GOES_date, ' ', (strmid(max_time(i),0,2) + ':' + strmid(max_time(i),2,2)), ' ', let_class(i), num_class(i)
endif
endfor ;i loop

;read data
readcol, 'calc/'+lbl+'_time.dat', fits_id, asc_time, nor_time, date1, format = 'f,f,a,a', /silent
readcol, +lbl+'_NL_length.dat',fits_id, NL_length, /silent
readcol, +lbl+'_ang_on_disc.dat', fits_id, ang_on_disc, /silent
readcol, +lbl+'_av_NL_grad.dat', fits_id,av_NL_grad, /silent
readcol, +lbl+'_eff_dist.dat', fits_id, eff_dist, /silent
readcol, +lbl+'_total_flux.dat', fits_id, total_flux, /silent
readcol, +lbl+'_WLSG.dat', fits_id, WLSG, /silent, format = 'd,f'

close, 2

;convert UT time to mjd
    ;mdi data time
  file_length2 = n_elements(date1) - 1
  MDI_mjd = dblarr(file_length2+1)
  for j = 0, file_length2 do begin
    file1 = date1(j) + ' ' + nor_time(j)
    jd_struc1 = anytim2jd(file1)
    MDI_mjd(j) = double((jd_struc1.int-2400000.5) +jd_struc1.frac)
  endfor
     ;GOES data
  readcol,  'calc/Flares.dat', rel_dt, rel_time, rel_flare, format = 'a,a,a', /silent
  file_length3 = n_elements(rel_dt)-1
  gtest = 0 & if file_length3 gt 0 then gtest = 1
  if gtest eq 1 then begin
    GOES_mjd = dblarr(file_length3+1)
    for k = 0, file_length3 do begin
      file2 = rel_dt(k) + ' ' + rel_time(k)
      jd_struc2 = anytim2jd(file2)
      print, jd_struc2
      jd_struc50 = anytim2jd(rel_dt(k) + ' ' + (strmid(rel_time, 0,2) + strmid(rel_time,2,2)))
      print, jd_struc50
      if double((jd_struc2.int-2400000.5) + jd_struc2.frac) le MDI_mjd(file_length2) then $
        GOES_mjd(k) = double((jd_struc2.int-2400000.5) + jd_struc2.frac)
    endfor
  endif

set_plot, 'ps'
device, /color

IF gtest EQ 1 THEN BEGIN
  GMJD = GOES_mjd(0)
  ;set up Verticle bars for flare location
  EDtop = 1.5*max(eff_dist) & EDbar = findgen(EDtop) & EDx = dblarr(n_elements(EDbar)) & EDx(*) = GMJD
  FLtop = max(total_flux) & FLbar = findgen(13) & FLx = dblarr(n_elements(FLbar)) & FLx(*) = GMJD
  GWtop = 1.5*max(WLSG) & GWbar = findgen(GWtop) & GWx = dblarr(n_elements(GWbar)) & GWx(*) = GMJD
  FLbar = FLbar*1E22
ENDIF

device, filename = 'plots/NL_length_vs_time.ps'
plot, MDI_mjd, NL_length, title = 'Neutral Line length', xtitle = 'TIME (mjd)', ytitle = '(Mm)'
if gtest eq 1 then oplot, GOES_mjd, NL_length(GOES_mjd), psym=2, symsize=2

device, filename = 'plots/ang_on_disc_vs_time.ps'
plot, MDI_mjd, ang_on_disc, title = 'Central Angle on Disc', xtitle = 'TIME (mjd)', ytitle = 'Degrees'
if gtest eq 1 then oplot, GOES_mjd, ang_on_disc, psym=2, symsize=2

device, filename = 'plots/av_NL_grad_vs_time.ps'
plot, MDI_mjd, av_NL_grad, title = 'Average Neutral Line Gradient'
if gtest eq 1 then oplot, GOES_mjd, av_NL_grad(GOES_mjd), psym=2, symsize=2

device, filename = 'plots/eff_dist_vs_time.ps'
plot, MDI_mjd, eff_dist, title = 'Effective Separation vs Time', xtitle = 'Date (Year 2000)', ytitle = '(pixels)', charthick=4,charsize=2, xthick=2, ythick=2, thick=3, xticks=5, xminor=9, xtickname=['3/16','3/17','3/18','3/19','3/20','3/21']
if gtest eq 1 then oplot, EDx, EDbar
oplot, [51621.981,51621.981], [0,10E24] ;flare 1
oplot, [51623.697,51623.697], [0,10E24] ;flare 2 


device, filename = 'plots/total_flux_vs_time.ps'
plot, MDI_mjd, total_flux, title = 'Total Unsigned Flux vs Time', xtitle = 'Date (Year 2000)', ytitle = '(Mx)', charthick=4,charsize=2, xthick=2, ythick=2, thick=3, xticks=5, xminor=9, xtickname=['3/16','3/17','3/18','3/19','3/20','3/21']
if gtest eq 1 then oplot, FLx, FLbar
oplot, [51621.981,51621.981], [0,10E24] ;flare 1
oplot, [51623.697,51623.697], [0,10E24] ;flare 2 


device, filename = 'plots/WLSG_vs_time.ps'
plot, MDI_mjd, WLSG, title = 'GWILL vs Time', xtitle = 'Date (Year 2000)', ytitle = '(Mm)', charthick=4,charsize=2, xthick=2, ythick=2,thick=3, xticks=5, xminor=9, xtickname=['3/16','3/17','3/18','3/19','3/20','3/21'] 
if gtest eq 1 then oplot, GWx, GWbar
oplot, [51621.981,51621.981], [0,10E24] ;flare 1
oplot, [51623.697,51623.697], [0,10E24] ;flare 2 

device,/close

end