;; GO Tutorial
;;----------------------------------------------------------

; Load in a file from the archive (data in ~SD/data/SDDB/):

archive,'sye',2004040218,2     ;-> radar ID, date+hour, duration (in hours)

;Loads the measured radar parameters into a 2-D array "range gate" x "beam"

;If you want more info on any proceedure or function just use: pri,'function', e.g.:

pri,'archive'

;-----------------------------------------------------------

print_file_info                               ;-> displays information on the currently loaded file
print_parameter_info                          ;-> displays information on the currently loaded parameters
list_parameters                               ;-> displays information on other valid parameters
set_parameter_list,'vel_err','pwr_s','elev'   ;-> specifies a different set of 3 parameters
archive,'sye',2004040218,2                    ;-> re-read the file loading in the new parameters
print_parameter_info                          ;-> displays information on the currently loaded parameters
default_parameter_list                        ;-> loads back the default list
archive,'sye',2004040218,2

set_scatter              ;-> dislays information on the different options regarding ground scatter

help,a1,a2,a3            ;-> a1,a2,a3 are the 3 parameter arrays read into from the FITACF file
help,a                   ;-> a is the "current" parameter used for plotting etc.

help,scan_id             ;-> scan_id is the CP_ID of the radar operating mode
print,scan_id(0)         ;-> e.g. 150 is common time data

help,tfreq               ;-> the transmit frequency of the radar (kHz)

help,beam_time           ;-> the time (in seconds since the start of the year) of the corresponding beam

; e.g. take a random beam, number 123:

print,time_convert('time',3,year=2004,sec=beam_time(123))

; gives the time in UT of the beam

;-----------------------------------------------------------

; Navigate to a time within the file

go_time,1814                 ;-> goes to 18:14:00 - the start time of scan #7

print,scan_mark(0:19)        ;-> labels the beams where each new scan starts
help,where(scan_mark eq 1)   ;-> corresponds to 60 scans in this file
print,where(beam_scan eq 7)  ;-> the beams corresponding to scan 7
print,beam_dir(123)          ;-> beam 123 in the file is beam direction 11 on scan #7

;------------------------------------------------------------

; Plot a polar projection

plot_polar                          ;-> plots a polar projection view of the radar data
overlay_polar_coast,/south          ;-> overlays the coastline

plot_polar,/clip                    ;-> set the map limits to fit the radar field of view
overlay_polar_beams                 ;-> plots the radar field of view
overlay_polar_beams,beams=[11,12]   ;-> plots beam direction 11

erase
plot_polar_grid                     ;-> plots a blank polar grid
overlay_polar_coast,/south
overlay_polar_beams,/fill,col=253   ;-> produce a plot showing just radar FOV
overlay_polar_beams

erase
plot_polar_grid                     
overlay_sd_fov,20051007,1800,hemi=-1  ;-> plot all radar FOVs

plot_polar_fov,20051007,/reimei         ;-> plot REIMEI footprints
plot_polar_fov,20051007,/camera         ;-> plot camera FOVs

;-----------------------------------------------------------

; Try other coordinate systems

set_coords                                  ;-> the current coordinate system

set_coords,'geog'                           ;-> set geographic coordinates
xmin=40 & xmax=120 & ymin=-80 & ymax=-50    ;-> set the map limits for plotting
plot_map                                    ;-> plot!

set_coords,'geog/stereo'
xmin=50 & xmax=90 & ymin=-80 & ymax=-50
plot_map

set_coords,'mag'
xmin=70 & xmax=120 & ymin=-85 & ymax=-65
plot_map

set_coords,'mag/stereo'
xmin=80 & xmax=120 & ymin=-90 & ymax=-65
plot_map

set_coords,'range'
plot_map

;-------------------------------------------------------------

; Display a range time intensity (RTI) plot

set_coords,'gate'    ;-> specify range gate as the y-axis
time,1800,2000       ;-> set the rti time limits
plot_rti

;-----------------------------------------------------------

;Write the figure to a postscript file

ps_start,'piccy.ps'      ;-> open the postscript file
set_format,/landscape    ;-> set the desired format
clear_page
plot_rti
ps_close                 ;-> close the postscript file
ps_preview               ;-> preview the postscript file

ps_start,'piccy.ps'
set_format,/portrait
clear_page
plot_rti
ps_close
ps_preview

;-----------------------------------------------------------

; Stack plot of time-series from individual range gates

clear_page
vel
map,y=[15,35]                        ;-> rti y-range

plot_rti_panel,1,2,0,0,beam=0       ;-> can specify multiple panels on a page
plot_title
plot_colour_bar,2,0

oplot,[18,20],[20,20],thick=2
oplot,[18,20],[24,24],thick=2
oplot,[18,20],[28,28],thick=2
plot_stack_panel,1,2,0,1,beam=0,range=[20,24,28],inc=500,/bar,y=[-500,1500]

;---------------------------------------------------------------

; Stereo data and different operating modes

auto_concat,'~adrian/example/DATA/','han',2007121006,2   ;-> to use fitacf files other than archived ones
pwr_l
set_scatter,0
set_coords,'range'
set_format,/landscape
plot_rti

print,beam_time(0:19)       ;-> can see that file contains multiple records from the same time
make_mono,153               ;-> to isolate one mode from stereo data
print,beam_time(0:19)
plot_rti

auto_concat,'~adrian/example/DATA/','han',2007121006,2
make_mono,-6401                                       ;-> select the other mode

plot_rti_panel,1,1,0,0
oplot,[7.005,7.005],[0,3000],linestyle=2,thick=3

;------------------------------------------------------------------

; Polar plots of different radar modes (i) myopic

auto_concat,'~adrian/example/DATA/','han',2007121006,2
make_mono,153
go_time,0658
clear_page
plot_polar_panel,2,1,0,0,/clip
plot_map_info,2,1,0,0,no_freq=0,/SQUARE
overlay_polar_coast
overlay_polar_beams

print,lagfr(0),smsep(0)             ;-> the transmitted pulse lag time to 1st range
                                    ;-> and range separation (in microseconds)
print,3e8*.001*lagfr(0)/2.*1e-6     ;-> in km
print,3e8*.001*smsep(0)/2.*1e-6     ;-> in km

plot_polar_panel,2,1,1,0,/clip,/no_data
auto_concat,'~adrian/example/DATA/','han',2007121006,2
make_mono,-6401
go_time,0702
overlay_polar_data
overlay_polar_coast
overlay_polar_beams
plot_map_info,2,1,1,0,no_freq=0,/SQUARE
plot_title
plot_colour_bar

print,lagfr(0),smsep(0)
print,3e8*.001*lagfr(0)/2.*1e-6
print,3e8*.001*smsep(0)/2.*1e-6

;-----------------------------------------------------------------

; Polar plots of different radar modes (ii) high-time-resolution

auto_concat,'~adrian/example/DATA/','han',2007120110,2
set_scatter,0
print,beam_dir(0:19)                     ;-> can see interleaved modes of scanning and camped beams

make_mono,153
goto_scan,7                              ;-> can navigate by scan number rather time
clear_page
plot_polar_panel,2,1,0,0,/clip
plot_map_info,2,1,0,0,no_freq=0,/SQUARE
overlay_polar_coast
overlay_polar_beams

auto_concat,'~adrian/example/DATA/','han',2007120110,2
make_mono,-26007
goto_scan,7
plot_polar_panel,2,1,1,0,/clip
plot_map_info,2,1,1,0,no_freq=0,/SQUARE
overlay_polar_coast
overlay_polar_beams
overlay_polar_beams,beams=[9,10]

plot_title
plot_colour_bar

;-------------------------------------------------------------------

; High time resolution RTI

clear_page

auto_concat,'~adrian/example/DATA/','han',2007120110,2
make_mono,153
print,beam_time((where(beam_dir eq 9))(1))-beam_time((where(beam_dir eq 9))(0))
;-> shows the time resolution of the data = 60 seconds

set_coords,'gate'
map,y=[30,50]
vel
set_scale,-200,200,10
plot_rti_panel,1,2,0,0,beam=9

auto_concat,'~adrian/example/DATA/','han',2007120110,2
make_mono,-26007
print,beam_time((where(beam_dir eq 9))(1))-beam_time((where(beam_dir eq 9))(0))
;-> shows the time resolution of the data = 3 seconds

set_coords,'gate'
map,y=[30,50]
set_scale,-200,200,10
plot_rti_panel,1,2,0,1,beam=9
plot_colour_bar
plot_title

;--------------------------------------------------------------------

; Large-scale convection

archive,'han',2001012114,2
go_time,1504
vel
set_scatter,2
plot_polar

;overlay data from many radars on 1 plot
;---------------------------------------

archive,'pyk',2001012114,2
go_time,1504 & overlay_polar_data

archive,'sto',2001012114,2
go_time,1504 & overlay_polar_data

archive,'gbr',2001012114,2
go_time,1504 & overlay_polar_data

archive,'kap',2001012114,2
go_time,1504 & overlay_polar_data

archive,'sas',2001012114,2
go_time,1504 & overlay_polar_data

archive,'kod',2001012114,2
go_time,1504 & overlay_polar_data

archive,'pgr',2001012114,2
go_time,1504 & overlay_polar_data

;--------------------------------------------------------------------

; Loading in the archived global convection pattern data

archive_mp,20010121,/north,forder=8   ;-> date, hemisphere, [fit order = 8, 10, 12 is optional, default = 8]
help,map_data,/str                    ;-> shows the contents of the map file structure
mappot_time,1504                      ;-> navigating a convection map file
overlay_potential_contours            ;-> overlay the potential contours on the current plot
fitv
losv                                  ;-> specify the type of vectors to plot
truev
overlay_vectors,/bw                   ;-> overlays the vectors on the plot [black/white optional]

plot_potential_data                   ;-> plotting just a basic convection map
plot_potential_data,/storm            ;-> lower the equatorward plot boundary

overlay_fov,'han',/noclip             ;-> overlay an individual radar field of view

;--------------------------------------------------------------------------------

; Spacecraft footprints

Re=6375.                      ;-> define earth radius in km
sc_pos=[-15*Re, 3*Re, 1*Re]   ;-> define the spacecraft position in GSE coordinates
inputs=[4.0,-20,-2.0,-7.0]    ;-> define the input parameters (see below)

;Input_param:
; pdyn - solar wind dynamic pressure in nPa (ideally between 0.5 and 10)
; dst - Dst index value (ideadlly between -100 and +20)
; BYIMF and BZIMF - IMF components in nT (ideally between -10 and +10)

;plot the mapped position
overlay_spacecraft,sc_pos,input_param=inputs,sc_size=2

;----------------------------------------------------------------------
;----------------------------------------------------------------------

; Convection mapping analysis

archive_mp,20000109,/north
mappot_time,2216
plot_potential_page,2,1   ;-> plots a page of x,y convection maps 

; plot individual convection map panels

clear_page
mappot_time,2216
plot_potential_panel,2,2,0,0
mappot_time,2218
plot_potential_panel,2,2,1,0

;----------------------------------------------------------------------

; Generate a map file from the archive of grid files
; can change:
; (1) The time limits of the analysis (archive is full day files)
; (2) The latitude of the convection boundary
; (3) The IMF components (and hence the statistical model)
; (4) The propagation delay used in adjusting for the ACE upstream location
; (5) The order of the spherical harmonic fit

gen_map_pot, 20000109, 2200, 2300, hmb=63, imfb=[2,3,4], ordr=6

mappot_time,2216
plot_potential_panel,2,2,0,1
mappot_time,2218
plot_potential_panel,2,2,1,1

; can see how to adjust the fitting parameters to improve the resulting pattern

;------------------------------------------------------------------------

clear_page
mappot_time,2216
plot_potential_panel,2,2,0,0
mappot_time,2218
plot_potential_panel,2,2,1,0

; Generate a grid file from the archive of fitacf files
gen_grid_file, 2000010922,2, int=1200, scanl=1200, chan='A'

;**note the long integration time**

; Then regenerate the map file from this new grid file
gen_map_pot, '20000109n.grid', 2200, 2300, hmb=63, ordr=6, imf_del='1:12'

;-> Notice this time that if we leave out the fixed IMF values,
;-> real data from the ACE s/c are read in 

mappot_time,2216
plot_potential_panel,2,2,0,1
mappot_time,2218
plot_potential_panel,2,2,1,1

; other parameters that can be set when generating map files are:
;;/south - looks for southern hemisphere files
;;ordr - set this to the order of fit you want, if omitted order=6
;;hmbnd - set this to fix the latitude of the Hepner-Maynard boundary at magnetic
;;        local midnight. If omitted a dynamic boundary is 
;;        determined using:
;;hm_vel - (default=100) and
;;hm_npnts - (default=3)
;;imf_del - a string 'hh:mm': set this to fix the IMF delay time, if omitted a variable 
;;  	    delay is determined from SW data
;;imfb - set this to a 3 element array to fix the IMF rather than using the ACE data
;;model_doping - set this to the doping level of the model convection patterns, where:
;;         0=light, 1=medium, 2=heavy, 3=extreme
;;verbose - set this keyword to print out detailed info about the fit
;;noread - set this keyword to inhibit loading in if the newly generated MAP file
;;

; other parameters that can be set when generating grid files are:
;;/south - looks for southern hemisphere files
;;int= - integrate the grid data into records of length "int" seconds
;;scanl= - ignore the scan flag and instead use a fixed length scan of "scanl" seconds
;;pmax= - exclude data with a power greater than pmax
;;pmin= - exclude data with a power less than pmin
;;vmax= - exclude data with an absolute velocity greater than vmax
;;vmin= - exclude data with an absolute velocity less than vmin
;;wmax= - exclude data with a spectral width greater than wmax
;;wmin= - exclude data with a spectral width less than wmin
;;vemax= - exclude data with a velocity error greater than vemax
;;vemin= - exclude data with a velocity error less than vemin
;;fmax= - set the maximum allowed variation in frequency to be fmax Hertz
;;chan= - process data from channel "A|B" for stereo mode data
;;minrng= - exclude data from range gates lower than minrange
;;maxrng= - exclude data from range gates higher than maxrange
;;/ion - exclude data marked as ground scatter (default)
;;/gs - exclude data marked as ionospheric scatter
;;/both - do not exclude data based on the type flag
;;/nav - do not perform temporal averaging
;;xtd - create extended output that includes both the power and spectral width
;;inertial - calculate the grid file in the inertial reference frame
;;ebm - exclude the beams listed in "ebm", a comma separated list of beam numbers

;-----------------------------------------------------------------------------------------------------------
;-----------------------------------------------------------------------------------------------------------

;STATISTICS

restore,'~SD/data/SDDB/ade/stat_files/stats_2001n.sav',/verb   ;-> restore the stats file
stats_n=stats
restore,'~SD/data/SDDB/ade/stat_files/stats_2001s.sav'
stats_s=stats
help,stats_s
help,stats_s(10000),/str

; NO ROUTINES EXIST FOR PLOTTING STATISTICS DATA
; JUST HAVE TO USE IDL TO DO WHAT YOU WANT!

; EXAMPLE:
; Show that the azimuthal sense of the convection over the polar cap
; is opposite in the two hemispheres

good_coverage=where(stats_s.n_points GE 400, count)   ;-> identify intervals of good data cov.
print,count

good_coverage=where(stats_s.n_points GE 400 AND $     ;-> identify simultaneous good intervals
  stats_n.n_points GE 400, count)                     ;-> in both hemispheres
print,count 

clear_page
plot,stats_n(good_coverage).ve(36),stats_s(good_coverage).ve(36),psym=1, $
  yrange=[-1000,1000],position=[0.2,0.05,0.8,0.95]
fit=poly_fit(stats_n(good_coverage).ve(36),stats_s(good_coverage).ve(36),1)
oplot,indgen(2000)-1000,fit(1)*(indgen(2000)-1000)+fit(0)

; etc. etc....