pyb_traj.py

"""
Functions used by pyBalloon to calculate balloon trajectories
"""

import numpy as np
import datetime

import pyb_aux
import pyb_gfs
import pyb_io

import param_file as p

#################################################################################################################

def update_files(used_weather_files=None, data=None, current_date=None, current_time=None, current_loc=None, total_time=[0], params=None, balloon=None, parachute=None):
	"""
	Update weather files to newest and best available (closest in time to current time)

	Arguments
	=========
	used_weather_files : dict
		Dictionary of weather files that have been used so far. The keys represent the times at which the weather files started being used.
	data : dictionary
		Dictionary containing weather forecast data plus properties calculated using the calc_properties() function
	current_datestr : string
		Current datestr in prediction
	current_time : float
		Current time in prediction
	current_loc : floats in tuple
		Current (latitude in degrees, longitude in degrees, altitude in km) in prediction
	total_time : list
		List of increments in time for steps in trajectory so far
	params : dictionary
	   	Dictionary of parameters determining how the trajectory is calculated, e.g. with interpolation, descent_only etc.
	balloon : dictionary
		Dictionary of balloon parameters, e.g. burtsradius, mass etc.
	parachute : dict
		Dictionary of prachute parameters, e.g. area, drag coefficient
	para
	Return:
		Data of new weather file(s), updated weather file keys, and updated used_weather_files.
	"""

	res = str(int(-4*params['resolution'] + 6))

	updated = False

	current_lat, current_lon, current_alt = current_loc

	weather_files = list(data.keys())
	weather_files.sort()
	keys = weather_files

	time_keys = [(int(int(file[15:19])/100.) + int(file[20:23])) for file in weather_files]
	time_hhh = [int(file[20:23]) for file in weather_files]

	max_datestr, max_time, max_hhhh, max_hhh = keys[-1][6:14], time_keys[-1], keys[-1][15:17], keys[-1][20:23]
	
	current_hours, current_minute, current_seconds = pyb_gfs.convert_hours(current_time)
	date_current = datetime.datetime(int(current_date[:4]), int(current_date[4:6]), int(current_date[6:]), \
					 current_hours, current_minute, current_seconds)

	sum_hour = int(max_hhhh)+int(max_hhh)
	rem_hours = sum_hour % 24
	extra_days = (sum_hour - rem_hours)/24.
	date_weather_file = datetime.datetime(int(max_datestr[:4]), int(max_datestr[4:6]), int(max_datestr[6:]), rem_hours, 0, 0)
	date_weather_file += datetime.timedelta(days=extra_days)

	data_lats, data_lons = data[keys[-1]]['lats'], data[keys[-1]]['lons']
	min_lat, max_lat, min_lon, max_lon = min(data_lats), max(data_lats), min(data_lons), max(data_lons)

	min_lat_diff, max_lat_diff = np.abs(min_lat-current_lat), np.abs(max_lat-current_lat)
	min_lon_diff, max_lon_diff = np.abs(min_lon-current_lon), np.abs(max_lon-current_lon)

	# edge of tile
	if current_lon-min_lon < 0 or max_lon-current_lon < 0:

		print('At %s, %.1f utc, passed the edge of tile, reading new tile from weather forecasts...' % (current_date, current_time))

		data[keys[-2]] = prepare_data(weather_file=keys[-2], loc0=current_loc, \
			current_time=current_time, balloon=balloon, params=params, parachute=parachute)

		data[keys[-1]] = prepare_data(weather_file=keys[-1], loc0=current_loc, \
			current_time=current_time, balloon=balloon, params=params, parachute=parachute)

		updated = True

		keys = list(data.keys())
		keys.sort()

	# past time of latest weather file
	elif date_current >= date_weather_file:

		print('Date and time is now %s, %.1f utc, adding new weather file...' % (current_date, current_time))

		new_weather_files = pyb_gfs.get_interpolation_gfs_files(datestr=current_date, utc_hour=current_time, params=params)	
		new_hhhh, new_hhh1 = int(int(new_weather_files[-1][15:19])/100), int(new_weather_files[-1][21:24])
		new_hhhh0, new_hhh01 = int(int(new_weather_files[0][15:19])/100), int(new_weather_files[0][21:24])

		# one weather file will be the same, but we want to use the new current location
		data[new_weather_files[-2]] = prepare_data(weather_file=new_weather_files[-2], loc0=current_loc, \
			current_time=current_time, balloon=balloon, params=params, parachute=parachute)

		data[new_weather_files[-1]] = prepare_data(weather_file=new_weather_files[-1], loc0=current_loc, \
			current_time=current_time, balloon=balloon, params=params, parachute=parachute)

		used_weather_files[current_time] = new_weather_files

		updated = True

		keys = list(data.keys())
		keys.sort()
		
	return data, keys, used_weather_files, updated

###############################################################################################################

def calc_time_frac(current_time=None, weather_files=None):
	"""
	Calculate how far the current time is from each interpolation file as fractions of difference in time over total time.

	E.g. if the current time is 10 then on the lhs we have 6 and on the rhs we have 12, 10 is then 2/3 of the way away from 6 and 1/3 away from 12.
	the fractions needed are then 1-2/3 for 6 and 1 - 1/3 for 12

	Arguments
	=========
	current_time : float
		Current time of trajectory
	weather_files : list
		List of the interpolation files currently being used
	Return:
		Two tuples containing the interpolation file name and corresponding fraction
	"""

	weather_files.sort()
	times = [(int(int(file[15:19])/100.) + int(file[20:23])) for file in weather_files]

	earlier_file = weather_files[-2]
	later_file = weather_files[-1]

	earlier_time = times[-2]
	later_time = times[-1]

	if later_time !=  24:
		dt1 = current_time - (earlier_time % 24)
		dt2 = (later_time % 24) - current_time
	else:
		dt1 = current_time - (earlier_time % 24)
		dt2 = later_time - current_time		

	dt_total = 3.0

	frac1 = 1. - dt1/dt_total
	frac2 = 1. - dt2/dt_total

	return (earlier_file, frac1), (later_file, frac2)

#################################################################################################################

def read_data(loc0=None, weather_file=None, params=None):
	"""
	Read in model forecast data

	Arguments
	=========
	loc0 : floats in tuple
		(latitude in degrees, longitude in degrees, altitude in km) of initial point
	weather_file : string
		Name of weather file from data is to be read
	params : dictionary
		Dictionary of parameters determining how the trajectory is calculated, e.g. with interpolation, descent_only etc.
	Return:
		Dictionary containing the forecast data
	"""

	lat0, lon0, alt0 = loc0

	in_dir = p.path + p.forecast_data_folder
	
	# note top, left, bottom, right ordering for area
	area = (lat0 + (params['tile_size']/2.), lon0 - (params['tile_size']/2.), lat0 - (params['tile_size']/2.), lon0 + (params['tile_size']/2.))

	print('Reading GFS data from ' + str(weather_file) + '.grb2...')

	file_dir = in_dir		
	model_data = pyb_io.read_gfs_single(directory=file_dir + weather_file, area=area, alt0=alt0, descent_only=params['descent_only'])[0]

	return model_data

#################################################################################################################

def movement2ll(lat_rad=None, lon_rad=None, alt=None, dx=None, dy=None):
	"""
	Calculate new lat/lon coordinates given original location and Cartesian displacements.

	Arguments
	=========
	lat_rad : float
		Current latitude in radians
	lon_rad : float
		Current longitude in radians
	alt : float
		Current altitude in meters
	dx : float
		East - west movement in meters (East is positive)
	dy : float
		North - south movement in meters (North is positive)
	Return:
		New coordinates: latitude [radians], longitude [radians], and distance traveled [km]
	"""

	# radius derived from WGS84 reference ellipsoid with altitude added to it.
	Re = pyb_aux.earth_radius(lat_rad)
	radius = Re + alt/1000. # convert alt to km
	
	# calculate distance travelled
	dist = np.sqrt(dx*dx + dy*dy)/1000. # Convert to km

	# calculate direction of movement
	theta = np.arctan2(dx, dy)

	cos_dr = np.cos(dist/radius)
	sin_dr = np.sin(dist/radius)
	sin_lat = np.sin(lat_rad)
	cos_lat = np.cos(lat_rad)
	cos_theta = np.cos(theta)
	sin_theta = np.sin(theta)

	# use haversine formula
	lat2 = np.arcsin(sin_lat * cos_dr + cos_lat * sin_dr * cos_theta)
	lon2 = lon_rad + np.arctan2(sin_theta * sin_dr * cos_lat, cos_dr - sin_lat * np.sin(lat2))

	return lat2, lon2, dist

#################################################################################################################

def calc_properties(data, weather_file, loc0, balloon, parachute, params):
	"""
	Calculate necessary properties, e.g. air densities and descent speeds, and carry the interpolation to more altitudes than in forecast data

	Arguments
	=========
	data : dictionary
		Dictionary containing data from forecast model read in by read_data()
	weather_file : string
		Name of weather file from which the data is to be read (if data is None)
	loc0 : floats in tuple
		(latitude in degrees, longitude in degrees, altitude in km) of initial point
	params : dictionary
		Dictionary of parameters determining how the trajectory is calculated, e.g. with interpolation, descent_only etc.
	balloon : dict
		Dictionary of balloon parameters, e.g. burtsradius, mass etc.
	parachute : dict
		Dictionary of prachute parameters, e.g. area, drag coefficient
	Return:
		Data appended with new properties
	"""

	lat0, lon0, alt0 = loc0

	data['air_densities'] = pyb_aux.air_density(data)

	if not params['descent_only']:

		data['balloon_radii'], gas_mass = pyb_aux.mooney_rivlin(data, balloon['radius_empty'], balloon['fill_radius'], balloon['thickness_empty'])
		data['balloon_volumes'] = pyb_aux.balloon_volume(data)
		total_mass = parachute['equip_mass'] + balloon['balloon_mass'] + gas_mass # kg
		data['lifts'] = pyb_aux.lift(data, total_mass)
		data['ascent_speeds'] = pyb_aux.ascent_speed(data, total_mass, balloon['Cd_balloon'])
		if 'simple_ascent_rate' in balloon:
			data['ascent_speeds'] = np.ones_like(data['ascent_speeds']) * balloon['simple_ascent_rate']

	data['descent_speeds'] = pyb_aux.descent_speed(data, parachute['equip_mass'], parachute['Cd_parachute'], parachute['parachute_areas'], params['altitude_step'], parachute['parachute_change_altitude'])

	data = pyb_aux.data_interpolation(data=data, alt0=alt0, step=params['altitude_step'], descent_only=params['descent_only'])

	return data

#################################################################################################################

def calc_displacements(data, balloon, params, parachute):
	"""
	Calculate the displacements in the x/y directions for all possible locations (lon/lat/alt) in weather data

	Arguments
	=========
	data : dictionary
		Dictionary containing weather forecast data plus properties calculated using the calc_properties() function
	params : dictionary
		Dictionary of parameters determining how the trajectory is calculated, e.g. with interpolation, descent_only etc.
	balloon : dict
		Dictionary of balloon parameters, e.g. burtsradius, mass etc.
	parachute : dict
		Dictionary of prachute parameters, e.g. area, drag coefficient
	Return:
		Data appended with additional calculated properties  		
	"""

 	# ascent properties		
	if not params['descent_only']:
		for i in range(len(data['ascent_speeds'])):
			data['ascent_speeds'][i] = np.array(list(0.5*data['ascent_speeds'][i][:-1] + 0.5*data['ascent_speeds'][i][1:]) + [data['ascent_speeds'][i][-1]])

	if not params['descent_only']:
		data['max_altitudes'], data['max_alt_idxs'] = pyb_aux.burst_altitude(data, balloon['burst_radius'])
		delta_t = params['altitude_step'] / data['ascent_speeds']
		data['ascent_time_steps'] = delta_t
		data['cumulative_ascent_times'] = np.cumsum(delta_t)/60.

 	# descent properties		
	for i in range(len(data['descent_speeds'])):
		data['descent_speeds'][i] = np.array(list(0.5*data['descent_speeds'][i][:-1] + 0.5*data['descent_speeds'][i][1:]) + [data['descent_speeds'][i][-1]])
		data['u_winds'][i] = np.array(list(0.5*data['u_winds'][i][:-1] + 0.5*data['u_winds'][i][1:]) + [data['u_winds'][i][-1]])
		data['v_winds'][i] = np.array(list(0.5*data['v_winds'][i][:-1] + 0.5*data['v_winds'][i][1:]) + [data['v_winds'][i][-1]])
	
	delta_t = -1*params['altitude_step'] / data['descent_speeds']

	data['descent_time_steps'] = np.array(delta_t)
	data['cumulative_descent_times'] = np.cumsum(delta_t)/60

	return data

#################################################################################################################

def calc_variable(grid_i, i, lon_rad, lat_rad, data, fracs, prop, resolution):
	"""
	Calculate the value of a variable at a given grid location/altitude

	Arguments
	=========
	grid_i : int
		index of list representing current lon/lat
	i : int
		index of list representing current altitude
	lon_rad : float
		Current longitude in radians
	lat_rad : float
		Current latitude in radians
	data : dictionary
		Dictionary containing weather forecast data plus properties calculated using the calc_properties() function
	fracs : list
		List containing the fractions used for the interpolation between weather files
	prop : string
		Name of variable to be evaluated
	resolution : float
		Resolution of the weather forecasts
	Return:
		Value of given variable  		
	"""

	t1, f1, t2, f2 = fracs
	keys = list(data.keys())
	data_lats, data_lons = np.radians(data[t1]['lats']), np.radians(data[t1]['lons'])

	x1, y1, x2, y2, low_left, up_left, low_right, up_right = pyb_aux.find_bilinear_points(grid_i=grid_i, i=i, lon_rad=lon_rad[-1], lat_rad=lat_rad[-1], \
						data_lons=data_lons, data_lats=data_lats, resolution=resolution)
	coords, inds = [x1, y1, x2, y2], [low_left, up_left, low_right, up_right]

	var = f1*pyb_aux.bilinear_interpolation(i=i, lon_rad=lon_rad[-1], lat_rad=lat_rad[-1], prop=data[t1][prop], coords=coords, inds=inds)
	var += f2*pyb_aux.bilinear_interpolation(i=i, lon_rad=lon_rad[-1], lat_rad=lat_rad[-1], prop=data[t2][prop], coords=coords, inds=inds)

	return var

#################################################################################################################

def calc_movements(data=None, used_weather_files=None, ini_conditions=None, params=None, balloon=None, parachute=None):
	"""
	Calculate the trajectory of the balloon/parachute given a start position & other input parameters

	Arguments
	=========
	used_weather_files : dictionary
		Dictionary of initial weather files used. The keys represent the times at which the weather files started being used (i.e. here the starting time)	
	ini_conditions : tuple of strings
		(Date of initial point, Initial time of trajectory, (latitude in degrees, longitude in degrees, altitude in km) of initial point
	params : dictionary
		Dictionary of parameters determining how the trajectory is calculated, e.g. with interpolation, descent_only etc.
	balloon : dict
		Dictionary of balloon parameters, e.g. burtsradius, mass etc.
	parachute : dict
		Dictionary of prachute parameters, e.g. area, drag coefficient
	Return:
		Output trajectory data, and dictionary of used_weather_files
	"""

	############################################################################################################

	# set general parameters and initial conditions
	datestr, utc_hour, loc0 = ini_conditions
	lat0, lon0, alt0 = loc0

	lat_rad, lon_rad, all_alts = [np.radians(lat0)], [np.radians(lon0)], [alt0]

	current_time = utc_hour
	current_loc = loc0
	
	keys = list(data.keys())
	time_key0 = (int(int(keys[0][15:19])/100.) + int(keys[0][20:23]))

	alts = data[keys[0]]['altitudes'] # alts, lats and lons are the same for all weather files (if we dont give it different areas)
	data_lats, data_lons  = np.radians(data[keys[0]]['lats']), np.radians(data[keys[0]]['lons'])
	size = int(np.sqrt(len(data_lats)))

	x_prop, y_prop = 'u_winds', 'v_winds'
	t_props = ['ascent_time_steps', 'descent_time_steps']
	speed_props = ['ascent_speeds', 'descent_speeds']

	tfutures, speeds = [], []
	total_time, dists, dists_u, dists_v = [0], [0], [0], [0]
	stage, props_index, i, max_i, timer = 1, 0, 0, 0, 0

	if params['drift_time'] == 0:
		stage_update = 2
	else:
		stage_update = 1

	if not params['descent_only']:
		print('Calculating ascent...')

	# calc trajectory
	while True:

		if not (stage == 1 and params['descent_only']):

			# determine fractions for the interpolation between forecasts
			(t1, f1), (t2, f2) = calc_time_frac(current_time=current_time, weather_files=keys)

			delta_t = current_time - (int(int(t1[15:19])/100.) + int(t1[20:23])) # difference between forecast and model time
			tfuture = f1*int(t1[20:23]) + f2*int(t2[20:23])

			# Find the closest grid point
			diff = np.sqrt((data_lats - lat_rad[-1])**2 + (data_lons - lon_rad[-1])**2)
			grid_i, = np.where(diff == diff.min())
			grid_i = grid_i[0]
			min_diff = diff[grid_i]

			dx = calc_variable(grid_i, i, lon_rad, lat_rad, data, (t1, f1, t2, f2), x_prop, params['resolution'])
			dy = calc_variable(grid_i, i, lon_rad, lat_rad, data, (t1, f1, t2, f2), y_prop, params['resolution'])

			if stage != 2:

				dt = calc_variable(grid_i, i, lon_rad, lat_rad, data, (t1, f1, t2, f2), t_props[props_index], params['resolution'])
				speed = calc_variable(grid_i, i, lon_rad, lat_rad, data, (t1, f1, t2, f2), speed_props[props_index], params['resolution'])

			else:

				dt = 5 # 60*60*6
				speed = 0
				
			dx *= dt
			dy *= dt

		if stage == 1:

			if not params['descent_only']:

				max_alt = f1*data[t1]['max_altitudes'][grid_i]+f2*data[t2]['max_altitudes'][grid_i]

				if all_alts[-1] >= max_alt:

					stage += stage_update
					props_index += 1

					final_i = 0

					print('Balloon burst at %.0f m' % (all_alts[-1]))

					continue

			else:

				# Mimick the movement during ascent if we only want descent
				if alts[i] >= alt0:
					stage += stage_update
					props_index += 1
					final_i = i
					continue

			i += 1
			max_i = i

		elif stage == 2:

			if timer == 0:
				print('Calculating drift trajectory...')

			timer += dt
			if timer >= params['drift_time']*60: # to seconds
				stage += 1
				continue

		elif stage == 3:

			if i == 0:

				(t1, f1), (t2, f2) = calc_time_frac(current_time=current_time, weather_files=keys)

				delta_t = current_time - (int(int(t1[15:19])/100.) + int(t1[20:23])) # difference between forecast and model time
				tfuture = f1*int(t1[20:23]) + f2*int(t2[20:23])
				tfutures.append(tfuture)

				print('Calculating descent...')

				break

			i -= 1

		if stage == 2:

			lat, lon, dist = movement2ll(lat_rad=lat_rad[-1], lon_rad=lon_rad[-1], alt=all_alts[-1], dx=dx, dy=dy)
			alt = all_alts[-1] # assume altitude does not change during drift

		elif stage != 2 and not (stage == 1 and params['descent_only']):

			lat, lon, dist = movement2ll(lat_rad=lat_rad[-1], lon_rad=lon_rad[-1], alt=alts[i], dx=dx, dy=dy)
			alt = alts[i]

		if not (params['descent_only'] and stage == 1):
			
			lon_deg = check_near_meridians(lon=np.degrees(lon), tile_size=params['tile_size'])
			current_loc = np.degrees(lat), lon_deg, alt
			lon = np.radians(lon_deg)

			speeds.append(speed)
			lat_rad.append(lat)
			lon_rad.append(lon)
			dists.append(dist)
			total_time.append(dt)
			all_alts.append(alt)
			tfutures.append(tfuture)

		# total_time[-1] is dt, i.e. always the same
		if current_time + total_time[-1]/3600. >= 24.:

			year, month, day = int(datestr[:4]), int(datestr[4:6]), int(datestr[6:])
			date = datetime.datetime(year, month, day) + datetime.timedelta(days=1)
			datestr = str(date.year) + str(date.month).zfill(2) + str(date.day).zfill(2)

		current_time = (float(utc_hour) + np.cumsum(np.array(total_time))[-1]/3600) % 24

		# update weather files
		data, keys, used_weather_files, updated = update_files(used_weather_files=used_weather_files, data=data, \
			current_date=datestr, current_time=current_time, current_loc=current_loc, \
			total_time=total_time, params=params, balloon=balloon, parachute=parachute)

		if updated:
			alts = data[keys[-1]]['altitudes'] # alts, lats and lons are the same for all weather files (if we dont give it different areas)
			data_lats, data_lons  = np.radians(data[keys[-1]]['lats']), np.radians(data[keys[-1]]['lons'])
			size = int(np.sqrt(len(data_lats)))

	speeds.append(calc_variable(grid_i, i, lon_rad, lat_rad, data, (t1, f1, t2, f2), speed_props[props_index], params['resolution']))

	try:
		# get more accurate end-point based on elevation data
		if params['check_elevation']:

			print('Getting new endpoint based on elevation...')

			elevation = pyb_aux.get_elevation(lon=np.degrees(lon_rad[-1]), lat=np.degrees(lat_rad[-1]))
			if np.abs(elevation - all_alts[-1]) > params['altitude_step']/10.:

				new_end_point, new_alt = pyb_aux.get_endpoint(data=(np.degrees(np.array(lat_rad)), np.degrees(np.array(lon_rad)), \
										    np.array(all_alts), np.array(dists)))

				diff = np.sqrt((np.degrees(np.array(lat_rad)) - new_end_point[0])**2 + (np.degrees(np.array(lon_rad)) - new_end_point[1])**2)
				index = np.where(diff == min(diff))[0][-1]

				lat_rad, lon_rad, all_alts = lat_rad[:index], lon_rad[:index], all_alts[:index], 
				dists, speeds, total_time = dists[:index+1], speeds[:index+1], total_time[:index+1]

				tfutures = tfutures[:index+1]

				lat_rad.append(np.radians(new_end_point[0]))
				lon_rad.append(np.radians(new_end_point[1]))
				all_alts.append(new_alt)
	except:
		print('Could not find endpoint based on elevation.\nPerhaps there is no data available for this region?')

	############################################################################################################

	# Convert the result array lists to Numpy 2D-arrays
	output = {}
	output['lats'] = np.degrees(np.array(lat_rad)) # to decimal degrees
	output['lons'] = np.degrees(np.array(lon_rad)) # to decimal degrees
	output['alts'] = np.array(all_alts)

	output['dists'] = np.array(dists)
	output['times'] = np.cumsum(np.array(total_time))/60 # to minutes
	output['distance'] = np.sum(np.array(dists))
	output['speeds'] = np.array(speeds)

	output['tfutures'] = np.array(tfutures)
	output['mean_direction'] = pyb_aux.calc_mean_travel_direction(lon0=lon0, lat0=lat0, end_lon=float(output['lons'][-1]), \
								      end_lat=float(output['lats'][-1]))

	print('\nFinished!\n')

	final_lat, final_lon, final_alt = output['lats'][-1], output['lons'][-1], output['alts'][-1]

	if final_lon > 180:
		final_lon -= 360
	if final_lon < -180:
		final_lon += 360

	if final_lon < 0:
		lon_add = 'W'
	else:
		lon_add = 'E'

	if final_lat < 0:
		lat_add = 'S'
	else:
		lat_add = 'N'

	# print out relevant quantities
	print('Maximum altitude: %.2f m' % np.max(all_alts))
	print('Landing location: lat=%.4f deg %s, lon=%.4f deg %s, alt=%.2f m' % (final_lat, lat_add, final_lon, lon_add, final_alt))
	print('Impact velocity: %.2f m/s' % output['speeds'][-1])
	print('Flight time: %.1f min' % (output['times'][-1]))
	print('Distance travelled: %.1f km' % output['distance'])
	print('Mean direction of travel: %.1f deg from North' % (90+(360-output['mean_direction'])))
	print('')

	return output, used_weather_files

#################################################################################################################

def prepare_data(weather_file=None, loc0=None, current_time=None, balloon=None, params=None, parachute=None):
	"""
	Prepare the data from a weather file for calculating the trajectory
	
	Arguments
	=========
	weather_file : string
		Name of weather file to be used and read
	loc0 : floats in tuple
		(latitude in degrees, longitude in degrees, altitude in km) of initial point
	current_time : float
		Current time of the trajectory	
	params : dictionary
		Dictionary of parameters determining how the trajectory is calculated, e.g. with interpolation, descent_only etc.
	balloon : dict
		Dictionary of balloon parameters, e.g. burtsradius, mass etc.
	parachute : dict
		Dictionary of prachute parameters, e.g. area, drag coefficient
	Return:
		Dictionary containing data ready for starting the trajectory calculations
	"""

	model_data1 = read_data(loc0=loc0, weather_file=weather_file, params=params)
	model_data2 = calc_properties(data=model_data1, weather_file=weather_file, loc0=loc0, balloon=balloon, params=params, parachute=parachute)
	model_data3 = calc_displacements(data=model_data2, balloon=balloon, params=params, parachute=parachute)

	return model_data3

#################################################################################################################

def check_near_meridians(lon, tile_size):
	"""
	If longitude near prime or anti meridian
	Make sure that we can grab a continuous tile from the gfs data

	Arguments
	=========
	lon : float
		Longitude to check
	tiles_size : float
		Size of tile we wish to read from the gfs file
	Returns:
		Longitude in 0 -> 360 range if too close to the dateline or
		longitude in -180 -> 180 range if too close to prime meridian
	"""

	# if close to dateline
	if lon < 0 and (lon - tile_size/2.) < -180:
		lon += 360

	# if close to prime meridian
	if lon > 0 and (lon + tile_size/2.) > 360:
		lon = (lon + 180) % 360 - 180

	return lon

#################################################################################################################

def run_traj(weather_files=None, ini_conditions=None, params=None, balloon=None, parachute=None):
	"""
	Run all functions to calculate the trajectory
	
	Arguments
	=========
	weather_files : list
		List of weather files used at the start
	ini_conditions : tuple of strings
		(Date of initial point, Initial time of trajectory, (latitude in degrees, longitude in degrees, altitude in km) of initial point
	params : dictionary
		Dictionary of parameters determining how the trajectory is calculated, e.g. with interpolation, descent_only etc.
	balloon : dict
		Dictionary of balloon parameters, e.g. burtsradius, mass etc.
	parachute : dict
		Dictionary of prachute parameters, e.g. area, drag coefficient
	Return:
		The calculated trajectories, list of fig dicionaries for each weather file, dictionary of used_weather_files, and list of time differences between trajectiry times and weather forecast files
	"""

	datestr, utc_hour, loc0 = ini_conditions
	lat0, lon0, alt0 = loc0

	loc0 = lat0, check_near_meridians(lon=lon0, tile_size=params['tile_size']), alt0
	ini_conditions = datestr, utc_hour, loc0

	data_array, used_weather_files = {}, {utc_hour : weather_files}
	for weather_file in weather_files:

		model_data = prepare_data(weather_file=weather_file, loc0=loc0, current_time=utc_hour, balloon=balloon, parachute=parachute, params=params)
		data_array[weather_file] = model_data
		del model_data

	trajectories, used_weather_files = calc_movements(data=data_array, used_weather_files=used_weather_files, ini_conditions=ini_conditions, \
		 balloon=balloon, parachute=parachute, params=params)

	return trajectories, used_weather_files

#################################################################################################################

def err_parallel(fit_params, hor_dist, tfut):
	"""
	Method to calculate error in direction of travel

	Arguments
	=========
	fit_params : tuple of floats
		Tuple of the parameters used to calculate the error
	hor_dist : float
		Total predicted horizontal distance travelled
	tfut : float
		Average t_future for the weather models used at each altitude step
	Return:
		Predicted error in direction of travel
	"""

	sigma_0, h, k, q = fit_params

	return np.sqrt(sigma_0**2 + h*hor_dist**2 + k*tfut**2)

#################################################################################################################

def err_perp(fit_params, hor_dist, tfut):
	"""
	Method to calculate error in direction perpendicular to direction of travel

	Arguments
	=========
	fit_params : tuple of floats
		Tuple of the parameters used to calculate the error
	hor_dist : float
		Total predicted horizontal distance travelled
	tfut : float
		Average t_future for the weather models used at each altitude step
	Return:
		Predicted error in direction perpendicular to direction of travel
	"""

	sigma_0, h, k, q = fit_params

	return np.sqrt(sigma_0**2 + h*hor_dist**2 + k*tfut**2)/q

#################################################################################################################

def determine_error(data=None):
	"""
	Method to predict the error on the end_point

	Arguments
	=========
	data : dictionary
		Dictionary containing trajectory data
	Return:
		Predicted errors in parallel and perpendicular direction wrt the direction of travel
	"""

	fit_params = (1.769442, 0.000316, 0.003567, 1.142628) # sigma_0, h, k, q: parameters determined from test flights

	tfut = np.mean(data['tfutures'])
	hor_dist = np.sum(np.array(data['dists']))
	parallel_err, perp_err = err_parallel(fit_params, hor_dist, tfut), err_perp(fit_params, hor_dist, tfut)

	print('The errors are: %.2f and %.2f km\n' % (parallel_err, perp_err))

	return parallel_err, perp_err

#################################################################################################################