pathGrid_devnotes.txt

pathGrid dev notes
Todd Anderson
24 May 2021

Contents:
(searchable tag)					(description)

1. pg_gridcross_gpu notes			development of pg_gridcross for GPU implementation
2. dispersion_analysis				extension of pathGrid project to analysis of sferic dispersions
	a. da_operation				 		plan for how this is going to work
	b. da_newstations					what can we gain from adding new WWLLN stations (e.g. at high latitudes)



==========================
pg_gridcross_gpu notes

May 24 Profile results (1 10-minute strokefile, on my laptop, using pg_gridcross_gpu inside pathGrid_10m
	pathGrid_10m:		31.160s		(100%)
	pg_gridcross_gpu:	31.079s		(99.7%)
	track2:				23.435s		(75.2%)

	-Is there any way to replace or reduce runtime of track2?
		-can reduce n_points in track2 from 400, but risk missing grid locations
		-on a cartesian grid, latitude (y) and longitude (x) vary monotonically over track; equivalent to rhumb line path
		-can we find points in an arbitrary coordinate system in which latitude and longitude vary monotonically between start and end points, then rotate those points into Earth coordinates?
			-e.g: rotate so start point (s1) is always at one pole (s1'); then the 400 points between start and end will vary only in latitude.
				-see rotation in spherical coordinates here: 
				https://stla.github.io/stlapblog/posts/RotationSphericalCoordinates.html
				-rotation axis defined by points on the equator at s1_longitude + 90 degrees and s1_longitude - 90 degrees
				-rotation angle = s1_colatitude = 90 degrees - s1_latitude
				-in frame where s1' is at a pole, grid locations traversed are 
					lats' = s1'_lat:grid_spacing:s2'_lat;
					lons' = s2'_lon*ones(length(lats'));
				-rotate back to Earth frame
			-how do we ensure that grid locations in Earth frame are captured in prime frame?
				


 
% strokelist: list of stroke-station paths
% 		columns:

%		time	strokelat	strokelon	stationlat	stationlon
% rows
% stroke-station1
% stroke-station2
% .
% .
% .
% stroke-stationN

grid_gpua = zeros(180,360,'gpuArray')

for j in N: 

	%generate 400-point track between stroke and station 	-->	track 		(2 x 400)

	%round points down to nearest integer lat/lon		 	--> track_round	(2 x 400)
	%get unique points only								 	-->	unique_tra 	(2 x U)
	%transform into indices 								--> unique_ind	(2 x U)
	%	[(-90,89),(-180,179)] --> [(1,180),(1,360)]


	% for each unique lat-lon point, increment
	for k in U:

		grid_gpua(unique_tra(1,k),unique_ind(2,k)) = grid_gpua(unique_tra(1,k),unique_ind(2,k)) + 1

	end

end  

=======================
dispersion_analysis

operation: how do we extend stroke-station path analysis to include dispersion?
- dispersion analysis allows us to parameterize ionosphere
	- i.e. based on Wait and Spies 2-parameter ionosphere model
	- note that Morris Cohen et al have a 4-parameter model for the D-region that uses lightning sferics for probe signals
- dispersion is calculated from Sfiles
	- each station generates an Sfile every (1? 10?) minutes
	- Sfiles contain waveform information for each sferic detected at that station within the file's time range
	- the matlab script SfilePlot.m (and function togafit.m) calculates dispersion from Sfiles
- procedure:
	0. File formats
		a. Sfile format:
			each line of Sfile is a sferic detected at that station, with following comma-separated data ("example"):
			Magic header 													("W210")
			wwlln site id 													("42")
			UTC time YYYY-MM-DDThh:mm:ss 									("2012-10-28T23:46:00")
			UTC toga offset in microseconds									("64171")
			RMS amplitude of waveforms [units???]							("482")
			TOGA offset in seconds from start of waveforms					("0.000423515")
			The three dispersion fit parameters fit0 , fit1 , << fit2;		("0.000122209,-22.342,995620")
			a dispersion fit ok flag (0 or 1)								("1")
			sampling frequency in Hz										("47999.46798")
			number of waveforms samples										("64")
			the waveform samples											("-0.00014624,0.000582488,-0.000494393, ...")
			--> Sfiles are generated at each station, once per minute
			--> name format: SyyyymmddHHMM, where time in filename specifies start time of file
		b. Rfile format:
			station ID | time after start of hour (seconds) | rms amplitude of e field at station (uncalibrated)
			example:
			10 			 60.085788 							 2709
			--> Rfiles are generated by the network once per 10 minutes
			--> name format: RyyyymmddHHMM, where time in filename specifies start time of file
			--> Rfile times are time of group arrival ("TOGA" or "toga")
			--> Rfile time + last hour time = Sfile UTC time (element 3) + UTC toga offset time (element 4)
			--> check: for a given station, lines in Rfiles and Sfiles with matching times should have identical RMS amplitudes
	1. for each stroke-station pair in APfile, find associated sferic in Sfile
		a. APfile.data (.mat version): each row is one lightning stroke with following columns:
			(1:6) 	Date: yyyy, mm, dd, HH, MM, SS.SSS
			(7) 	latitude: fractional degrees north
			(8)		longitude: fractional degrees east
			(9)		residual fit error in mircoseconds
			(0)		number of WWLLN stations contributing to location
		b. APfile.power (.mat version): each column is one lightning stroke. Column N in APfile.power corresponds to Row N in APfile.data:
			(odd rows)	station ID (see stations.dat for list of station IDs)
			(even rows) energy of stroke detected by station listed in previous row, in arbitrary sound card units
			--> it is convenient to replace all even row values with NaN, -1, or other flag that cannot be mistaken for station ID
		c. stations.dat, or stations.mat, lists station IDs, station names, latitude, longitude
		d. for each station with Sfiles in analysis time range:
			1. find all instances of that station ID in AP.power (remember to exclude even rows!)
				--> index = Nx2 vector of [row, col] pairs corresponding to each instance of station ID
			2. for each instance of station ID, find corresponding row in AP.data
				--> data_stID = AP.data(index(:,2))
			3. for each row in data_stID, model propagation time from stroke to station
				- station location given in stations.mat
				--> time_station = time_stroke + distance(lat_stroke, lon_stroke, lat_station, lon_station, 'units')./c_vlf;
					where	time_station is time of sferic being recorded at station in question
							time_stroke is time of lightning stroke (i.e. time given in APfile)
							distance() is a MATLAB function that estimates distance
								lat_stroke, lon_stroke is the stroke location
								lat_station, lon_station is the station location from stations.mat
								'units' is the output units (e.g. km)
								may be other arguments!
							c_vlf is the propagation speed of VLF radio waves in the Earth-ionosphere waveguide, approximately equal to the speed of light
		e. get list of relevant sferic arrival times from Sfile
			QUESTION: do we want TOGA, or start of sferic?  probably TOGA
			1. for each row in Sfile.data, get time in format comparable to APfile time format
				--> time = timetransform(Sfile.data(:,3))		where timetransform() is a function that converts Sfile time format to desired time format
			2. add TOGA offsets to sferic times
				--> toga = Sfile.time + Sfile.data(:,4)
		f. match times between APfile.data_stID and Sfile.toga
			1. for each row in APfile.data_stID, find min(APfile.data_stID - Sfile.toga)
				--> need both index and value
				--> if abs(min(...)) > threshold, no match
				--> else, have matched each line APfile containing certain stID to a toga at that station
	2. for each stroke-station pair, calculate sferic dispersion measured by station
		a. 
	3. map stroke-station paths a la pathGrid.grid_crossings
		note this distribution will include only paths ending at stations for which Sfiles are available
	4. plots
		a. non-binned stroke-station paths colored by dispersion
		b. lat-lon binned stroke-station paths colored by:
			1. mean dispersion
			2. difference in mean dispersion from last hour average
	5. calculate 2-parameter ionosphere
	6. plots
		a. ionosphere parameter maps
	7. 


new stations: what can we gain from adding new stations (or new kinds of stations) to WWLLN?
- arctic stations
	try plotting s-s paths with simulated arctic stations
	what factors are important for number of new stations + station placement?
		distance to existing lightning stations and lightning sources
		ice sheets -- sferics don't propagate well over Greenland IS

other data sources
- ionosondes
- radars e.g. EISCAT, EISCAT_3D
	 not sure if data are publicly available, but you can request observation time