;+
; NAME:
;	forecast
; PURPOSE:
;	Updates all the plots and synoptic maps for the IPS forecast Web pages
;	using the output files from the last run of the tomography program.
; CATEGORY:
;	pro/tool
; CALLING SEQUENCE:
	PRO forecast, utnow, $
		raw_files		= raw_files		, $
		nocorrelation	= nocorrelation	, $
		fast			= fast			, $
		filter			= filter		, $
		magnetic		= magnetic		, $
		gif				= gif			, $
		png				= png			, $
		env_filter		= env_filter	, $
		nosave			= nosave		, $
		ftp				= ftp			, $
		silent			= silent		, $
		_extra			= _extra
; INPUTS:
;	utnow			current time as time structure or Carrington variable
;					if not defined, it is set to the system time
; OPTIONAL INPUT PARAMETERS:
;	/raw_files		process 'raw' tomography files, i.e. use files directly
;					output by the tomography files in the $NAGOYA/slow/raw
;					and $NAGOYA/fast/raw directories.
;					By default, the processed files in the slow/final and
;					fast/final directories are used.
;	/nocorrelation	skips production of the fore- and aftcast correlation
;					graphs. These take a long time to produce for the
;					time-dependent tomography, and are usually created only
;					once a day (when the tomography is run)
;	/fast			processes ime-dependent tomography files
;					(by default the corotating tomography files are processed)
;	filter=filter	scalar; type: string
;						either 'nv3d' or 'nv3f';
;						overrides prefix value set by href=vu_type_insitu=
;	/magnetic		process magnetic field, instead of plasma
;						density and velocity
;	/png, /gif		sets image format for output images (png is the default)
;	/nosave			suppresses output of graphics as GIF files
;						(used for debugging)
;	env_filter=env_filter
;					NOT SURE ANYMORE WHY THIS IS DONE
;					scalar; type: string; default: none
;					name of an env var to which the value of 'filter' is
;					assigned.
;	/ftp			OBSOLETE
;					forces html files and images to be copied to a remote
;					webserver (used to be cassfos02)
; OUTPUTS:
;	(GIF files; HTML files)
; INCLUDE:
	@compile_opt.pro							; On error, return to caller
; CALLS:
;	InitVar, IsType, Carrington, TimeSet
;	TimeGet, TimeUnit, TimeOp, TimeSystem, vu_localskymap
;	vu_synopticmap, destroyvar, forecast_info, forecast_ftp
;	vu_earthskymap, vu_remoteview, hide_env, vu_insitu
;	vu_select, vu_type_insitu, vu_type_skymap, vu_get
;	vu_gettime, vu_prefix, forecast_cfg
; PROCEDURE:
;	The tomography data should be located in $NAGOYAslow. Output files
;	(GIF files and updated HTML files) are written to $NAGOYA/slow/image
; MODIFICATION HISTORY:
;	MAR-2000, Paul Hick (UCSD/CASS)
;	JUL-2002, Paul Hick (UCSD/CASS)
;		Suppressed time axis on synoptic map when time-dependent tomography
;		files are processed.
;		Fixed a couple of problems related to processing of output from time-
;		dependent tomography. The identification of the time of the last
;		tomography run was wrong. The correlation plots now use e3 instead of
;		ea files.
;	SEP-2002, Paul Hick (UCSD/CASS)
;		Adapted to accept magnetic field 'bb3d' files
;		Added /forcecd keyword to vu_select
;	NOV-2002, Paul Hick (UCSD/CASS)
;		Added keyword /nowrap to one call of vu_gettime (calculating the
;		3D matrices for corotating tomography displays)
;	APR-2003, Paul Hick (UCSD/CASS)
;		Completed switch to png files (gif still available through keyword /gif)
;	JUL-2003, Paul Hick (UCSD/CASS)
;		Removed calculation of correlations from ea* and e3* files (this produced
;		the cast_insitu_****_v_c_n_c.gif images)
;	AUG-2007, Paul Hick (UCSD/CASS)
;		Introduced Carrington calls (replacing arg_time)
;	APR-2013, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Reorganized verbosity of all messages.
;		Removed /forcecd keyword everywhere (FindAllSubDirs now ignores this)
;-
InitVar, silent		, 0

InitVar, fast		, /key
InitVar, ftp		, /key
InitVar, raw_files	, /key
InitVar, nosave		, /key
InitVar, magnetic	, /key
InitVar, nocorrelation	, /key

InitVar, gif		, /key
InitVar, png		, /key
png OR= 1-gif

speed   = 1-magnetic
density = 1-magnetic
brad    =   magnetic
btang   =   magnetic


kind = (['slow','fast'])[fast]

dump	= 1-nosave
local2ut= forecast_cfg(/local2ut)
maxelo  = forecast_cfg(/maxelo, /degrees)
upto    = forecast_cfg(/round_time)
rounds  = round(forecast_cfg(/round_time, /mmsec))/1000
scale   = forecast_cfg(/scale_skymap)
fmt_rot	= forecast_cfg(/crformat)

delt	= forecast_cfg(/insitu_delt)
trange	= forecast_cfg(/insitu_trange)
tstep	= forecast_cfg(/insitu_tstep)

volsize	= [0.0, 1.5] ;forecast_cfg(/remote_volsize)
ngrid   = forecast_cfg(/remote_ngrid)
fovsize = forecast_cfg(/remote_fovsize)
color_v	= [0, 100] ;forecast_cfg(/remote_color_v)
color_n	= [125, 225] ;forecast_cfg(/remote_color_n)
range_v	= [200.0, 400.0] ; forecast_cfg(kind, /remote_range_v)
range_n	= [4.0, 20.0]    ;forecast_cfg(kind, /remote_range_n)
opa_v	= [0.3, 0.0]     ;forecast_cfg(kind, /remote_opacity_v)
opa_n	= [0.0, 0.5]     ;forecast_cfg(kind, /remote_opacity_n)
power_v	= 3.0 ; forecast_cfg(kind, /remote_opacity_power_v)
power_n	= 2.0 ; forecast_cfg(kind, /remote_opacity_power_n)
planar_breakn  = forecast_cfg(/planarcut_break_n)
planar_breakv  = forecast_cfg(/planarcut_break_v)

ylshow  = forecast_cfg(kind, /ylshow, /degrees)

; If set to 'fast' or 'slow' then sources overplotted on skysweep maps
; are extracted from a file nagoya.fast.good or nagoya.slow.good
; instead of from the nagoya.yyyy files. This file is created by the
; tomography program, and is moved to $NAGOYA/daily (see sync_ips_daily).

use_good_sources = fast ? 'fast' : ''

; raw_dir: directory where raw tomography files are located
; dat_dir: directory where tomography files for images are stored
; img_dir: ouput dir where gifs and html files are written by IDL
; sav_dir: output dir for backup of gif files (the images are put in subdirs)
; htm_dir: WWW dir. Gifs and html files are copied from img_dir to htm_dir.

nagoya	= forecast_cfg(/root)
root	= filepath(root=nagoya, subdir=kind, '')

raw_dir   = filepath(root=root, subdir='raw'  , '')
final_dir = filepath(root=root, subdir='final', '')
dat_dir = raw_files ? raw_dir : final_dir

ice_dir = dat_dir	;magnetic ? dat_dir : filepath(root=dat_dir, subdir='nice', '')


img_dir = filepath(root=root, subdir='image', '')
sav_dir = filepath(root=root, subdir='image', '')
htm_dir = filepath(root=forecast_cfg(/html), subdir=kind,'current')

; If destination does not exist then output is to the terminal

IF dump THEN destination = [img_dir, sav_dir]

uday = TimeUnit(/day)
umin = TimeUnit(/minute)
usec = TimeUnit(/sec)

; Get forecast time from input, or use system time.

live_run = IsType(utnow,/undefined)

CASE live_run OF
0: utnow = Carrington(utnow,/get_time)
1: BEGIN
	utnow = TimeSystem(local2ut,/silent)
	utnow = TimeGet(utnow, usec, roundt=rounds)		; Round to nearest hour
	say, threshold=2, silent=silent, $
		'now is '+TimeGet(utnow, /ymd, upto=umin)+	$
		' (CR'+string(Carrington(utnow),format=fmt_rot)+')'
END
ENDCASE

location = forecast_cfg(utnow, /remote_location)

; Check when last email from Nagoya arrived
;mail = filepath(root=getenv('DAT'), subdir=['nagoya','daily'], 'logmail.ips')

; Type of data (nv or bb) to be read

InitVar, filter, vu_prefix(magnetic=magnetic,/silent)
IF strpos(filter,'*') NE -1 THEN BEGIN
	say, threshold=2, silent=silent, $
		[	"file filter with wildcard specified, '"+filter+"'", $
			'... ABORTING'	]
	RETURN
ENDIF

IF IsType(env_filter,/defined) THEN BEGIN
	setenv, env_filter+'='+filter
	add2label = filter
ENDIF

prefix = filter
filter += '*.*'

; First find the time of the last tomography run preceding utnow by checking
; the files present in raw_dir. This time is needed (even if the final files
; are used to construct the images) for the construction of the skysweep
; maps and to update the html explanations.

; For the corotating tomography we could extract the times from the file names,
; but not for the time-dependent tomography. So we read the headers.
; The vu_select call sorts entries into chronological order. For time-dependent
; tomography all files for the same times (from different runs) are also sorted
; on marker value.

marker   = fast
highmark = fast

hdrs = vu_select(	$
		path	= raw_dir	, $	; Check raw data
		filter	= filter	, $	; File filter (nv3d*.*, etc.)
		marker	= marker	, $	; If set, only get files with markers
		highmark= highmark	, $	; If set, only get files with highest marker
		rots	= times		, $	; UT time (NOT forecast time)
		rounds	= 1800		, $	; Rounds 'times' array
		imark	= imark		, $	; marker array
		count	= count		, $	; # files
		silent	= silent	)

IF count EQ 0 THEN BEGIN
   say, threshold=2, silent=silent, $
		'no raw files '+hide_env(filepath(root=raw_dir, filter))
   RETURN
ENDIF

files = vu_get(hdrs, /file) 							; Extract file names from headers
forecast_times = vu_get(hdrs, /forecast_time)			; Extract forecast times from headers

destroyvar, hdrs

; For corotating tomography 'last_run' refers to the last run.
; For time-dependent tomography it points to the last time step of
; the last run. The entire time-dependent run is identified
; by the marker value 'imark' (and the same forecast time).

last_run = n_elements(forecast_times)-1

last_forecast_time = forecast_times[last_run]			; Time of last run
IF IsType(imark, /defined) THEN imark = imark[last_run]	; Only for time-dependent tomography

IF live_run THEN	$
	say, threshold=2, silent=silent,  $
		[	'for '+TimeGet(last_forecast_time, /ymd, upto=umin)		+ $
			' ('+string(Carrington(last_forecast_time),format=fmt_rot)+')'	+ $
			' ~'+strmid(files[last_run],strlen(nagoya))	]

; Now go get the file(s) from dat_dir needed to construct all images.

CASE raw_files OF

0: BEGIN

	; The final data are averages based on the raw tomography output
	; calculated at regular time intervals. For the corotating tomography
	; the files are currently (May 2013) created at 6-hour intervals.
	; For the time-dependent tomography they exactly match the time steps
	; output by the tomography itself: currently (July 2002) at 6-hour
	; intervals.

	; Get the times from the file names; do not read the headers yet.
	; The array contains Carrington variables for a regular sequence
	; of times at the whole hour.

	destroyvar, files
	times = vu_select(files, path=dat_dir, filter=filter, /nohdr, count=count, rounds=1800, silent=silent)

	IF count NE 0 AND magnetic THEN BEGIN

		times_nv = vu_select(files_nv, path=dat_dir, filter='nv3*', /nohdr, rounds=1800, silent=silent)

		tmp = where_common( TimeGet(TimeOp(/subtract,times   ,times[0]),/diff,/full,uday),	$
							TimeGet(TimeOp(/subtract,times_nv,times[0]),/diff,/full,uday), inref=inref, count=count)
		IF count EQ 0 THEN BEGIN
			say, threshold=2, silent=silent, $
				[	'no common plasma and magnetic data files'	, $
					'... ABORTING'	]
			RETURN
		ENDIF
		destroyvar, times_nv

		; Only use the common files

		times = times[tmp]
		files = files[tmp]
		files_nv = files_nv[inref]

	ENDIF

	; Select all files within 15 days of utnow. We need these files
	; to construct the insitu times series for the time-dependent tomography

	IF count NE 0 THEN BEGIN
		tmp = TimeGet(TimeOp(/subtract,times,utnow),/diff,/full,uday)
		tmp = where( -15 LE tmp AND tmp LE 15, count)
	ENDIF

	IF count EQ 0 THEN BEGIN
		say, threshold=2, silent=silent, $
			[	'no final files within 15 days of '+TimeGet(utnow,/ymd)	, $
				'... ABORTING'	]
		RETURN
	ENDIF

	times = times[tmp]
	files = files[tmp]
	IF magnetic THEN files_nv = files_nv[tmp]

END

1: BEGIN

	; If we are working from the raw data (/raw_files set) we already have the
	; necessary information: we already have the last tomography run.
	; For corotating tomography only a single file will be left.

	IF NOT fast THEN BEGIN
		tmp = where( TimeGet(TimeOp(/subtract,forecast_times,last_forecast_time),/diff,/full,uday) EQ 0, count) ; Collect the last run
		IF count GT 1 THEN BEGIN
			say, threshold=2, silent=silent, $
				[	'multiple corotating files with same forecast time'	, $
					'... ABORTING'	]
			RETURN
		ENDIF
		tmp = tmp[0]
		times = times[tmp]					; File names and times for the last run
		files = files[tmp]
	ENDIF

	IF magnetic THEN BEGIN

		; When processing magnetic field files, require that there are a matching
		; number of plasma files available.

		tmp = fast ? 'nv3h*.*_'+string(format='(I5.5)',imark) : 'nv3d'+strmid(GetFileSpec(files,from='name'),strlen(prefix))
		files_nv = file_search( filepath(root=raw_dir,tmp+'.gz'),count=count_nv )

		IF count_nv NE count THEN BEGIN
			say, threshold=2, silent=silent, $
				[	'# magnetic files do not match # of plasma files'	, $
					'    # magnetic files is '+strcompress(count   ,/rem)	, $
					'    # plasma   files is '+strcompress(count_nv,/rem)	, $
					'... ABORTING'	]
			RETURN
		ENDIF

	ENDIF

END

ENDCASE

;===========
; KLUDGE:

; Offset the region of interest (shifting westward).
; This should reduce the amount of bad data in the shifted ROI somewhat.

roi_offset = forecast_cfg(/roi_offset)*(1-fast)

; END KLUDGE
;===========


tmp = vu_gettime(times, ut=utnow, bracket=bracket)

IF NOT fast THEN BEGIN
	IF n_elements(bracket) EQ 1 THEN bracket = bracket[0]

	; For corotating tomography find the two files that bracket utnow.
	; Only read these two files (this is not strictly necessary; it just cuts
	; down on the number of files to be read by vu_select).

	times = times[bracket]
	files = files[bracket]
	IF magnetic THEN files_nv = files_nv[bracket]

	bracket = n_elements(bracket) EQ 1 ? 0 : [0,1]
ENDIF

hdr = vu_select(files, count=count, /read, /get_roi, ff=ff, roi_offset=roi_offset, silent=silent+1, /check)

IF magnetic THEN BEGIN

	; The following arrays are only used for the remote views for magnetic data
	; (showing the plasma data with the current sheet embedded).

	; The remote views require data at time utnow only. Extract the magnetic data
	; bracketing utnow from hdr and ff, and read the matching plasma data.
	; Then interpolate to time utnow.

	hdr_bb = hdr[        bracket]
	ff_bb  = ff [*,*,*,*,bracket]
	hdr_nv = vu_select(files_nv[bracket], /read, /get_roi, ff=ff_nv, roi_offset=roi_offset, silent=silent+1, /check)

	; Note that keyword /nowrap is NOT included.

 	tmp    = ff_bb
	hdr_bb = vu_gettime(hdr_bb, tmp, ut=utnow, ff_ut=ff_bb)
	;hdr_bb = vu_gettime(hdr_nv, tmp, ut=utnow, ff_ut=ff_bb)

 	tmp    = ff_nv
	hdr_nv = vu_gettime(hdr_nv, tmp, ut=utnow, ff_ut=ff_nv)

ENDIF

IF NOT fast THEN BEGIN

	; For corotating tomography find the matrix at time utnow by interpolation.
	; Note that keyword /nowrap is present.

	tmp = ff
	hdr = vu_gettime(hdr, tmp, ut=utnow, count=count, ff_ut=ff, /nowrap)
ENDIF

IF count EQ 0 THEN RETURN						; No data read

; Make the plots

say, threshold=2, silent=silent, '-- Synoptic maps'

type = vu_type_insitu('synoptic', _extra=_extra, add2label=add2label,	$
	speed=speed, density=density, brad=brad, btang=btang)

vu_synopticmap, hdr, ff, $
	ut0			= utnow			, $
	type		= type			, $
	destination	= destination	, $
	png			= png			, $
	gif			= gif			, $
	upto		= upto			, $
	/earth						, $
	/forecast					, $
	nmajor		= 5				, $
	timeaxis	= 1-fast		, $
	silent		= silent		, $
	ylshow		= ylshow		, $
	/degrees					, $
	/strictlabels ;, from=uday

info = type[0].display+'_'+kind+type[0].label+type[1].label
forecast_info, kind, img_dir, info, utnow=utnow, time=last_forecast_time, upto=upto, silent=silent
info = type[0].display+'_'+kind
forecast_info, kind, img_dir, info, utnow=utnow, time=last_forecast_time, upto=upto, silent=silent+1

say, threshold=2, silent=silent, '-- Solar disk maps'

type = vu_type_insitu('earth_solardisk', _extra=_extra, add2label=add2label,	$
	speed=speed, density=density, brad=brad, btang=btang)

vu_solardisk, hdr, ff, $
	ut0			= utnow			, $
	type		= type			, $
	destination	= destination	, $
	png			= png			, $
	gif			= gif			, $
	upto		= upto			, $
	silent		= silent

info = type[0].display+'_'+kind+type[0].label+type[1].label
forecast_info, kind, img_dir, info, utnow=utnow, time=last_forecast_time, upto=upto, silent=silent
info = type[0].display+'_'+kind
forecast_info, kind, img_dir, info, utnow=utnow, time=last_forecast_time, upto=upto, silent=silent+1

IF magnetic THEN BEGIN
	say, threshold=2, silent=silent, '-- Synoptic Map at 0.07 AU (Solar Surface)'

	type = vu_type_insitu('synoptic', _extra=_extra, add2label=add2label,   $
		brad=brad)
	type[0].label = type[0].label + '_ss'

	vu_synopticmap, hdr, ff, $
		ut0			= utnow			, $
		type		= type			, $
		destination	= destination	, $
		png			= png			, $
		gif			= gif			, $
		upto		= upto			, $
		radius		= 0.07			, $
		/forecast					, $
		nmajor		= 5				, $
		timeaxis	= 1-fast		, $
		silent		= silent		, $
		ylshow		= ylshow		, $
		/degrees					, $
		/strictlabels ;, from=uday
ENDIF

say, threshold=2, silent=silent, '-- Ecliptic cross section'

type = vu_type_insitu('planarcut', _extra=_extra, add2label=add2label,	$
	speed=speed, density=density, brad=brad, btang=btang)
if speed   then type[0].breakval = planar_breakv
if density then type[1].breakval = planar_breakn

body = [jpl_body(/earth,/string)+'&', big_body('Stereo A'), big_body('Stereo B')]; ,jpl_body(/mars,/string)]+'&'

vu_planarcut, hdr, ff, $
	ut0			= utnow			, $
	type		= type			, $
	destination	= destination	, $
	png			= png			, $
	gif			= gif			, $
	upto		= upto			, $
	silent		= silent		, $
	/strictlabel				, $
	/from_centered_ecliptic		, $
	/heliocentric				, $
	body		= body			, $
	radius		= 1.5

info = type[0].display+'_'+kind+type[0].label+type[1].label
forecast_info, kind, img_dir, info, utnow=utnow, time=last_forecast_time, upto=upto, silent=silent
info = type[0].display+'_'+kind
forecast_info, kind, img_dir, info, utnow=utnow, time=last_forecast_time, upto=upto, silent=silent+1

; Time series at Earth in comparison with ACE data

say, threshold=2, silent=silent, '-- Time series @ Earth'

tm_label =	(['corotating','time-dependent'])[fast]+	$
			(['',' ('+strupcase(strmid(filter,0,3))+'/NOAA)'])[magnetic]

type = vu_type_insitu('earth_insitu', _extra=_extra, add2label=add2label,	$
	speed=speed, density=density, brad=brad, btang=btang)

source  = (['acesw','aceb'])[magnetic]

vu_insitu, hdr, ff, $
	ut0			= utnow				, $
	type		= type[0]			, $
	/forecast						, $
	destination	= destination		, $
	png			= png				, $
	gif			= gif				, $
	source		= source			, $
	/nofuture						, $
	/timeseries						, $
	delt		= delt				, $
	ts_range	= trange			, $
	ts_step		= tstep				, $
	chars		= type[0].charsize	, $
	thick		= type[0].thick		, $
	upto		= upto				, $
	tm_label	= tm_label			, $
	silent		= silent			, $
	digits		= 3*speed+3*brad	, $
	from		= uday

source  = (['rcelias','aceb'])[magnetic]

vu_insitu, hdr, ff, $
	ut0			= utnow				, $
	type		= type[1]			, $
	/forecast						, $
	destination	= destination		, $
	png			= png				, $
	gif			= gif				, $
	source		= source			, $
	/nofuture						, $
	/timeseries						, $
	delt		= delt				, $
	ts_range	= trange			, $
	ts_step		= tstep				, $
	chars		= type[1].charsize	, $
	thick		= type[1].thick		, $
	upto		= upto				, $
	tm_label	= tm_label			, $
	silent		= silent			, $
	digits		= 3*density+3*btang	, $
	from		= uday

info = type[0].display+'_'+kind+type[0].label+type[1].label
forecast_info, kind, img_dir, info, utnow=utnow, time=last_forecast_time, upto=upto, silent=silent
info = type[0].display+'_'+kind
forecast_info, kind, img_dir, info, utnow=utnow, time=last_forecast_time, upto=upto, silent=silent+1

; Time series at planets
; To make sure these use the same input data as the
; Earth time series, keep this block of code here!)

body = [	$
		jpl_body(/mercury,/string), $
		jpl_body(/venus  ,/string), $
		jpl_body(/mars	 ,/string), $
		;big_body('Ulysses'	 	 ), $
		big_body('Stereo A'	 	 ), $
		big_body('Stereo B'		 )  $
	]

FOR nbody=0,n_elements(body)-1 DO BEGIN

	say, threshold=2, silent=silent, '-- Time series @ '+body[nbody]

	tm_label = (['corotating','time-dependent'])[fast]
	destroyvar, source		; vu_insitu decides on source

	string_id = strlowcase( body[nbody] )+'_insitu'

	type = vu_type_insitu(string_id, _extra=_extra, add2label=add2label,	$
		speed=speed, density=density, brad=brad, btang=btang)

	vu_insitu, hdr, ff, $
		ut0			= utnow				, $
		type		= type[0]			, $
		/forecast						, $
		destination	= destination		, $
		png			= png				, $
		gif			= gif				, $
		source		= source			, $
		/nofuture						, $
		/timeseries						, $
		delt		= delt				, $
		ts_range	= trange			, $
		ts_step		= tstep				, $
		chars		= type[0].charsize	, $
		thick		= type[0].thick		, $
		upto		= upto				, $
		tm_label	= tm_label			, $
		silent		= silent			, $
		digits		= 3*speed+3*brad	, $
		body		= body[nbody]		, $
		from		= uday

	vu_insitu, hdr, ff, $
		ut0			= utnow				, $
		type		= type[1]			, $
		/forecast						, $
		destination	= destination		, $
		png			= png				, $
		gif			= gif				, $
		source		= source			, $
		/nofuture						, $
		/timeseries						, $
		delt		= delt				, $
		ts_range	= trange			, $
		ts_step		= tstep				, $
		chars		= type[1].charsize	, $
		thick		= type[1].thick		, $
		upto		= upto				, $
		tm_label	= tm_label			, $
		silent		= silent			, $
		digits		= 3*density+3*btang	, $
		body		= body[nbody]		, $
		from		= uday

	info = type[0].display+'_'+kind+type[0].label+type[1].label
	forecast_info, kind, img_dir, info, utnow=utnow, time=last_forecast_time, upto=upto, silent=silent
	info = type[0].display+'_'+kind
	forecast_info, kind, img_dir, info, utnow=utnow, time=last_forecast_time, upto=upto, silent=silent+1

ENDFOR

; Skymaps and skysweeps are skipped for magnetic field data

IF NOT magnetic THEN BEGIN

	say, threshold=2, silent=silent, '-- Sky snapshots'

	type1 = vu_type_skymap('earth_skysnap', /ips_speed )
	type2 = vu_type_skymap('earth_skysnap', /ips_glevel)

	vu_earthskymap, hdr, ff, $
		ut0			= utnow			, $
		type		= type1			, $
		destination	= destination	, $
		png			= png			, $
		gif			= gif			, $
		upto		= upto			, $
		maxelo		= -1			, $
		silent		= silent		, $
		/legend						, $
		/strictlabels				, $
		/compass					, $
		minelo		= 11.5			, $
		/degrees

	vu_earthskymap, hdr, ff, $
		ut0			= utnow			, $
		type		= type1			, $
		destination	= destination	, $
		png			= png			, $
		gif			= gif			, $
		upto		= upto			, $
		maxelo		= maxelo		, $
		scale		= scale			, $	; Passed to plotearthskymap
		silent		= silent		, $
		/legend						, $
		/strictlabels				, $
		/compass					, $
		minelo		= 11.5			, $
		/degrees

	vu_earthskymap, hdr, ff, $
		ut0			= utnow			, $
		type		= type2			, $
		destination	= destination	, $
		png			= png			, $
		gif			= gif			, $
		upto		= upto			, $
		maxelo		= -1			, $
		silent		= silent		, $
		/legend						, $
		/strictlabels				, $
		/compass					, $
		minelo		= 11.5			, $
		/degrees

	vu_earthskymap, hdr, ff, $
		ut0			= utnow			, $
		type		= type2			, $
		destination	= destination	, $
		png			= png			, $
		gif			= gif			, $
		upto		= upto			, $
		maxelo		= maxelo		, $
		scale		= scale			, $ ; Passed to plotearthskymap
		silent		= silent		, $
		/legend						, $
		/strictlabels				, $
		/compass					, $
		minelo		= 11.5			, $
		/degrees

	info = type1.display+'_'+kind+type1.label+type2.label
	forecast_info, kind, img_dir, info, utnow=utnow, time=last_forecast_time, upto=upto, silent=silent
	info = type1.display+'_'+kind
	forecast_info, kind, img_dir, info, utnow=utnow, time=last_forecast_time, upto=upto, silent=silent+1

	say, threshold=2, silent=silent, '-- Sky sweeps'

	type1 = vu_type_skymap('earth_skysweep', /ips_speed)
	type2 = vu_type_skymap('earth_skysweep', /ips_glevel)

	vu_localskymap, hdr, ff, $
		ut0			= last_forecast_time		, $
		type		= type1			, $
		destination	= destination	, $
		png			= png			, $
		gif			= gif			, $
		upto		= upto			, $
		silent		= silent		, $
		/legend						, $
		/compass					, $
		minelo		= 11.5			, $
		/degrees					, $
		/strictlabels				, $
		use_good_sources = use_good_sources	, $
		utnoon		= utnoon			; Output

	vu_localskymap, hdr, ff, $
		ut0			= last_forecast_time		, $
		type		= type1			, $
		destination	= destination	, $
		png			= png			, $
		gif			= gif			, $
		upto		= upto			, $
		silent		= silent		, $
		/legend						, $
		/compass					, $
		scale		= scale			, $ ; Passed to plotearthskymap
		maxelo		= maxelo		, $
		minelo		= 11.5			, $
		/degrees					, $
		/strictlabels				, $
		use_good_sources = use_good_sources

	vu_localskymap, hdr, ff, $
		ut0			= last_forecast_time		, $
		type		= type2			, $
		destination	= destination	, $
		png			= png			, $
		gif			= gif			, $
		upto		= upto			, $
		silent		= silent		, $
		/legend						, $
		/compass					, $
		minelo		= 11.5			, $
		/degrees					, $
		/strictlabels				, $
		use_good_sources = use_good_sources

	vu_localskymap, hdr, ff, $
		ut0			= last_forecast_time		, $
		type		= type2			, $
		destination	= destination	, $
		png			= png			, $
		gif			= gif			, $
		upto		= upto			, $
		silent		= silent		, $
		/legend						, $
		/compass					, $
		scale		= scale			, $ ; Passed to plotearthskymap
		maxelo		= maxelo		, $
		minelo		= 11.5			, $
		/degrees					, $
		/strictlabels				, $
		use_good_sources = use_good_sources

	; Argument utnoon is output from vu_localskymap,
	; so the forecast_info call must stay here

	info = type1.display+'_'+kind+type1.label+type2.label
	forecast_info, kind, img_dir, info, utnow=utnow, time=last_forecast_time, utnoon=utnoon, upto=upto, silent=silent
	info = type1.display+'_'+kind
	forecast_info, kind, img_dir, info, utnow=utnow, time=last_forecast_time, utnoon=utnoon, upto=upto, silent=silent+1

ENDIF

; For magnetic data the hcs is added to the density/velocity remote views.

say, threshold=2, silent=silent, '-- Remote view from behind Earth'

type = vu_type_insitu('earth_remoteview',/speed, /density, _extra=_extra)

CASE magnetic OF
0: BEGIN

	vu_remoteview, hdr, ff, $
		ut0				= utnow				, $
		type_ff			= type[0]			, $
		destination		= destination		, $
		png				= png				, $
		gif				= gif				, $
		view_time		= utnow				, $
		location		= location			, $
		/ecliptic							, $
		ngrid			= ngrid				, $
		volsize			= volsize			, $
		full_range		= range_v			, $
		full_opa		= opa_v				, $
		charsize		= type[0].charsize	, $
		upto			= upto				, $
		silent			= silent			, $
		voxel			= 2					, $
		opacity_power	= power_v			, $
		cue_radial		= [1,1.45,0.1,1]
		;fovsize=fovsize, /eqangle, full_color=color_v

	vu_remoteview, hdr, ff, $
		ut0				= utnow				, $
		type_ff			= type[1]			, $
		destination		= destination		, $
		png				= png				, $
		gif				= gif				, $
		view_time		= utnow				, $
		location		= location			, $
		/ecliptic							, $
		ngrid			= ngrid				, $
		volsize			= volsize			, $
		full_color		= color_n			, $
		full_range		= range_n			, $
		full_opa		= opa_n				, $
		charsize		= type[1].charsize	, $
		upto			= upto				, $
		silent			= silent			, $
		voxel			= 2					, $
		opacity_power	= power_n			, $
		cue_radial		= [1,1.45,0.1,1]
		;fovsize=fovsize, /eqangle

	info = type[0].display+'_'+kind+type[0].label+type[1].label
	forecast_info, kind, img_dir, info, utnow=utnow, time=last_forecast_time, upto=upto, silent=silent

	say, threshold=2, silent=silent, '-- Remote view from ecliptic north'

	; The '&' symbol adds orbits.
	body = [jpl_body(/mercury,/string), $
			jpl_body(/venus  ,/string), $
			jpl_body(/mars   ,/string)] + '&' ;, $
			;big_body('Ulysses')]+'&'
	; add stereo, but not the orbits.
	body = [body, big_body('Stereo A'), big_body('Stereo B')]

	vu_remoteview, hdr, ff, $
		ut0				= utnow				, $
		type_ff			= type[1]			, $
		destination		= destination[0]	, $
		png				= png				, $
		gif				= gif				, $
		view_time		= utnow				, $
		location		= [BadValue(1.0d0),!dpi/2,3.0d0]	, $
		/ecliptic							, $
		ngrid			= ngrid				, $
		/eqangle							, $
		fovsize			= fovsize			, $
		volsize			= volsize			, $
		full_color		= color_n			, $
		full_range		= range_n			, $
		full_opa		= opa_n				, $
		voxel			= 2					, $
		body			= body				, $
		xysize			= 275				, $
		module			= 'polar_splash'	, $
		silent			= silent	;, thick=2

	body = [body, jpl_body(/earth, /string)+'&']
	type = vu_type_insitu('planarcut', _extra=_extra, add2label=add2label,	$
		density=density)

	if density then type[0].breakval = planar_breakn
	;if density then type[0].breakval = gridgen(10,range=[2,10])

	vu_planarcut, hdr, ff, $
		ut0			= utnow				, $
		type		= type				, $
		destination	= destination[0]	, $
		png			= png				, $
		gif			= gif				, $
		xysize		= [300,300]			, $
		upto		= upto				, $
		silent		= silent			, $
		module		= 'ecliptic_cut'	, $
		/strictlabel					, $
		/from_centered_ecliptic			, $
		/heliocentric					, $
		body		= body				, $
		radius		= 2.0

END
1: BEGIN

	; terrible hack to adjust ff_bb to match the dimensions of ff

	nvdims = size(ff_nv,/dim)
	bbdims = size(ff_bb,/dim)
	IF total(nvdims EQ bbdims) NE n_elements(nvdims) THEN BEGIN
		hdr_bb.nlng = hdr_nv.nlng
		hdr_bb.nlat = hdr_nv.nlat
		hdr_bb.nrad = hdr_nv.nrad
		hdr_bb.distance_step = hdr_nv.distance_step
		ff_bi = fltarr(nvdims[0],nvdims[1],nvdims[2],bbdims[3])
		FOR i=0,2 DO $
			ff_bi[*,*,*,i] = interpolate(ff_bb[*,*,*,i], /grid, $
				(bbdims[0]-1)*gridgen(nvdims[0],/one), $
				(bbdims[1]-1)*gridgen(nvdims[1],/one), $
				(bbdims[2]-1)*gridgen(nvdims[2],/one))
		ff_bb = ff_bi
		; fingers crossed
	ENDIF

	type_bb = vu_type_insitu('earth_remoteview', /brad, add2label=add2label, _extra=_extra)
	vu_remoteview, hdr_nv, ff_nv, $
		ut0			= utnow				, $
		type_ff		= type[0]			, $
		destination	= destination		, $
		png			= png				, $
		gif			= gif				, $
		view_time	= utnow				, $
		location	= location			, $
		/ecliptic						, $
		ngrid		= ngrid				, $
		/eqangle						, $
		fovsize		= fovsize			, $
		volsize		= volsize			, $
		full_color	= color_v			, $
		full_range	= range_v			, $
		full_opa	= opa_v				, $
		charsize	= type[0].charsize	, $
		/meshhcs						, $
		hdr_bb		= hdr_bb			, $
		bb			= ff_bb				, $
		type_bb		= type_bb			, $
		upto		= upto				, $
		silent		= silent			, $
		voxel		= 2

	type_bb = vu_type_insitu('earth_remoteview', /btang, add2label=add2label, _extra=_extra)

	vu_remoteview, hdr_nv, ff_nv, $
		ut0			= utnow				, $
		type_ff		= type[1]			, $
		destination	= destination		, $
		png			= png				, $
		gif			= gif				, $
		view_time	= utnow				, $
		location	= location			, $
		/ecliptic						, $
		ngrid		= ngrid				, $
		/eqangle						, $
		fovsize		= fovsize			, $
		volsize		= volsize			, $
		full_color	= color_n			, $
		full_range	= range_n			, $
		full_opa	= opa_n				, $
		charsize	= type[1].charsize	, $
		/meshhcs						, $
		hdr_bb		= hdr_bb			, $
		bb			= ff_bb				, $
		type_bb		= type_bb			, $
		upto		= upto				, $
		silent		= silent			, $
		voxel		= 2

	type = vu_type_insitu('earth_remoteview',/brad, /btang, add2label=add2label)

	info = type[0].display+'_'+kind+type[0].label+type[1].label
	forecast_info, kind, img_dir, info, utnow=utnow, time=last_forecast_time, upto=upto, silent=silent
	info = type[0].display+'_'+kind
	forecast_info, kind, img_dir, info, utnow=utnow, time=last_forecast_time, upto=upto, silent=silent+1

END
ENDCASE

; One day forecast at Earth; need to read new data, so destroy hdr and ff

IF NOT nocorrelation THEN BEGIN

	destroyvar, hdr, ff

	say, threshold=2, silent=silent, '-- One-day forecast correlation'

	type = vu_type_insitu('cast_correlation', _extra=_extra, add2label=add2label,	$
		speed=speed, density=density, brad=brad, btang=btang)

	source = (['acesw','aceb'])[magnetic]

	; Changing raw_dir to ice_dir to read in more files for correlation

	vu_insitu_raw, hdr, ff, $
		ut0			= utnow				, $
		type		= type[0]			, $
		destination	= destination		, $
		png			= png				, $
		gif			= gif				, $
		source		= source			, $
		/nofuture						, $
		/correlation					, $
		delt		= delt				, $
		chars		= type[0].charsize	, $
		thick		= type[0].thick		, $
		upto		= upto				, $
		silent		= silent			, $
		cast_offset	= 1.0d0				, $
		ntrail		= 30				, $			; Use last 30 markers
		/strictlabels					, $
		path		= ([dat_dir,ice_dir])[fast], $
		filter		= filter

	source = (['rcelias','aceb'])[magnetic]

	vu_insitu_raw, hdr, ff, $
		ut0			= utnow				, $
		type		= type[1]			, $
		destination	= destination		, $
		png			= png				, $
		gif			= gif				, $
		source		= source			, $
		/nofuture						, $
		/correlation					, $
		delt		= delt				, $
		chars		= type[1].charsize	, $
		thick		= type[1].thick		, $
		upto		= upto				, $
		silent		= silent			, $
		cast_offset	= 1.0d0				, $			; Use last 30 markers
		ntrail		= 30				, $
		/strictlabels					, $
		path		= ([dat_dir,ice_dir])[fast], $
		filter		= filter

	; Five-day aftcast at Earth; need to read new data
	; For the corotating tomography fore- and aftcast are based on the same
	; files so read new data only for time-dependent tomography.

	say, threshold=2, silent=silent, '-- Five-day aftcast correlation'

	IF fast THEN destroyvar, hdr, ff

	type = vu_type_insitu('cast_correlation', _extra=_extra, add2label=add2label,	$
		speed=speed, density=density, brad=brad, btang=btang )

	source = (['acesw','aceb'])[magnetic]

	vu_insitu_raw, hdr, ff, $
		ut0			= utnow				, $
		type		= type[0]			, $
 		destination	= destination		, $
		png			= png				, $
		gif			= gif				, $
		source		 =source			, $
		/nofuture						, $
		/correlation					, $
		delt		= delt				, $
 		chars		= type[0].charsize	, $
		thick		= type[0].thick		, $
		upto		= upto				, $
 		silent		= silent			, $
		cast_offset	= -5.0d0			, $
		ntrail		= 30				, $			; Use last 30 markers
 		path		= ([dat_dir,ice_dir])[fast], $
		filter		= filter

	source = (['rcelias','aceb'])[magnetic]

 	vu_insitu_raw, hdr, ff, $
		ut0			= utnow				, $
		type		= type[1]			, $
 		destination	= destination		, $
		png			= png				, $
		gif			= gif				, $
		source		= source			, $
		/nofuture						, $
		/correlation					, $
		delt		= delt				, $
 		chars		= type[1].charsize	, $
		thick		= type[1].thick		, $
		upto		= upto				, $
		silent		= silent			, $
		cast_offset	= -5.0d0			, $
		ntrail		= 30				, $			; Use last 30 markers
		path		= ([dat_dir,ice_dir])[fast], $
		filter		= filter

	info = type[0].display+'_'+kind+type[0].label+'_c'+type[1].label+'_c'
	forecast_info, kind, img_dir, info, utnow=utnow, time=last_forecast_time, upto=upto, silent=silent
	info = type[0].display+'_'+kind
	forecast_info, kind, img_dir, info, utnow=utnow, time=last_forecast_time, upto=upto, silent=silent+1

	; Plot comparison of raw timeseries and spacecraft data.
	; We need to read raw data files, so destroy hdr and ff

	destroyvar, hdr, ff

	say, threshold=2, silent=silent, '-- Raw vs insitu correlation'

	type = vu_type_insitu('raw_correlation', _extra=_extra, add2label=add2label,	$
		speed=speed, density=density, brad=brad, btang=btang)

	source = (['acesw','aceb'])[magnetic]

	vu_insitu_raw, hdr, ff, $
		ut0			= utnow				, $
		type		= type[0]			, $
		/forecast						, $
		destination	= destination		, $
		png			= png				, $
		gif			= gif				, $
		source		= source			, $
		/nofuture						, $
		/correlation					, $
		delt		= delt				, $
		chars		= type[0].charsize	, $
		thick		= type[0].thick		, $
		upto		= upto				, $
		silent		= silent			, $
		/last_raw						, $
		/strictlabels					, $	
		path		= ([dat_dir,ice_dir])[fast], $
		filter		= filter

	source = (['rcelias','aceb'])[magnetic]

	vu_insitu_raw, hdr, ff, $
		ut0			= utnow				, $
		type		= type[1]			, $
		/forecast						, $
		destination	= destination		, $
		png			= png				, $
		gif			= gif				, $
		source		= source			, $
		/nofuture						, $
		/correlation					, $
		delt		= delt				, $
		chars		= type[1].charsize	, $
		thick		= type[1].thick		, $
		upto		= upto				, $
		silent		= silent			, $
		/last_raw						, $
		/strictlabels					, $
		path		= ([dat_dir,ice_dir])[fast], $
		filter		= filter

	say, threshold=2, silent=silent, '-- Raw vs insitu time series'

	type = vu_type_insitu('raw_timeseries', _extra=_extra, add2label=add2label, $
		speed=speed, density=density, brad=brad, btang=btang)

	source = (['acesw','aceb'])[magnetic]

	vu_insitu_raw, hdr, ff, $
		ut0			= utnow				, $
		type		= type[0]			, $
		/forecast						, $
		destination	= destination		, $
		png			= png				, $
		gif			= gif				, $
		source		= source			, $
		/nofuture						, $
		/timeseries						, $
		delt		= delt				, $
		chars		= type[0].charsize	, $
		thick		= type[0].thick		, $
		upto		= upto				, $
		silent		= silent			, $
		/last_raw						, $
		/strictlabels					, $
		path		= ([dat_dir,ice_dir])[fast], $
		filter		= filter

	source = (['rcelias','aceb'])[magnetic]

	vu_insitu_raw, hdr, ff, $
		ut0			= utnow				, $
		type		= type[1]			, $
		/forecast						, $
		destination	= destination		, $
		png			= png				, $
		gif			= gif				, $
		source		= source			, $
		/nofuture						, $
		/timeseries						, $
		delt		= delt				, $
		chars		= type[1].charsize	, $
		thick		= type[1].thick		, $
		upto		= upto				, $
		silent		= silent			, $
		/last_raw						, $
		/strictlabels					, $
		path		= ([dat_dir,ice_dir])[fast], $
		filter		= filter

	info = type[0].display+'_'+kind+type[0].label+type[1].label
	forecast_info, kind, img_dir, info, utnow=utnow, time=last_forecast_time, upto=upto, silent=silent
	info = type[0].display+'_'+kind
	forecast_info, kind, img_dir, info, utnow=utnow, time=last_forecast_time, upto=upto, silent=silent+1

ENDIF

; Copy images (png or gif) and html files to html directory
; Ftp to cassfos02 is requested (/ftp set)

IF getenv('DRYRUN') NE '1' THEN BEGIN

	img_type = '*.'+(['png','gif'])[0*png+1*gif]

	say, threshold=0, silent=silent, 'copy files from '+img_dir+' to '+htm_dir

	spawn, 'cp '+filepath(root=img_dir,'*.*')+' '+htm_dir
	IF ftp THEN forecast_ftp, kind, filepath(root=img_dir,img_type),	$
		filepath(root=img_dir,'*.html'), silent=silent

ENDIF ELSE BEGIN

	say, threshold=0, silent=silent, 'DRYRUN: not copying files from '+img_dir+' to '+htm_dir

ENDELSE

RETURN  &  END