Calculating distance travelled along a flight trajectory

Macro CalculateDistances()

Silent(1)
PauseUpdate

Variable i
Make/O/N=(numpnts(LAT)-1) dist_calc; dist_calc = NaN
do
	dist_calc[i] = (sin((LAT[i]-LAT[i+1])*pi/180/2))^2 + cos(LAT[i]*pi/180)*cos(LAT[i+1]*pi/180)*(sin((LON[i]-LON[i+1])*pi/180/2))^2
	dist_calc[i] = 2*atan2(sqrt(dist_calc[i]),sqrt(1-dist_calc[i]))
	i += 1
while(i < numpnts(LAT)-1)

//Save as both meters and feet.
Duplicate/O dist_calc dist_m
dist_m = dist_calc * 6371*1000	//convert to meters
dist_calc *= 3959*5280			//convert to feet

//Integrate along the path, such that the last point in dist_m_INT will be the total length travelled in meters.
Duplicate/O dist_m dist_m_INT
Integrate dist_m/D=dist_m_INT

End

Variable LAT_ORIGIN = [value]

Variable LON_ORIGIN = [value]

i+1

i

astrostu wrote: Does this help? It's code I wrote that I'm pretty sure uses Vincenty's formulae. The LAT is the latitude wave, LON is longitude wave. It calculates the distances between each successive point, storing them to "dist_calc" and "dist_m". It's not elegant, but it does the job.

*****IGOR CODE****
If you wanted to modify this to do it from a fixed point each time, then you would set a

Variable LAT_ORIGIN = [value]

and

Variable LON_ORIGIN = [value]

, change the LAT[i] and LON[i] to those origin variables, change

i+1

to

i

, and run it. You'd want to get rid of the do-while loop in that case, using Igor's matrix notation, but again, without that it'd still get the job done.

Hi

While searching for a solution to this myself I came across this thread. I'll add my own solutions here for completeness if anyone has need for it. 

Firstly, a function that takes lat,lon, plus optional start and end points, and altitude. The total distance can be calculated in 3 ways:

1. total distance over WGS-84 ellipsoid

2. total distance over WGS-84 ellipsoid accounting for altitude using simple pythagorean theorem - beware this works OK if you have high time resolution positional data, but not necessarily if you have two points far apart.

3. total spherical distance

// Calculate total distance along a set of lat-lon co-ordinates
// method=0 (default) - total distance over WGS-84 ellipsoid accounting for altitude above if given
// method=1 - total spherical distance 
Function DistanceAlongTrack(lat,lon[,stPnt,endPnt,method,alt])
	Wave lat,lon
	Variable stPnt,endPnt,method
	Wave alt
	
	// if no altitude is given make a zero wave in its place
	if (paramIsDefault(alt))
		Make/O/FREE/N=(dimSize(lat,0)) alt=0
		endif
	
	// Lat and Lon waves should be the same size 
	if (dimSize(lat,0)!=dimSize(lon,0))
		return 0
		endif
	
	// if no endpoint is given 	
	if (paramisDefault(endPnt))
		endPnt=dimSize(lat,0)
		endif

	Variable i,totDist
	for (i=stPnt;i<(endPnt-1);i+=1)
		if (method==0) // total distance over WGS-84 ellipsoid accounting for altitude above
			totDist+=EllipsoidalDistanceBetweenTwoPoints(lat[i],lon[i],lat[i+1],lon[i+1],alt1=alt[i],alt2=alt[i+1])
		elseif (method==1) // total spherical distance 		
			totDist+=SphericalDistanceBetweenTwoPoints(lat[i],lon[i],lat[i+1],lon[i+1])
		endif
	endfor
	
	// distance in meters	
	return totDist
end

 Secondly there is code to calculate total spherical distance using the Haversine formula

// calculates spherical distance between two points on earth given coords in decimal degrees
// uses haversine formula
// radius of earth is taken to be 6371 km (IUGG/WGS-84 mean)
Function SphericalDistanceBetweenTwoPoints(lat1, long1, lat2, long2)
    Variable lat1, long1, lat2, long2
    
   	// convert inputs in degrees to radians
	lat1*= pi/180
	long1*= pi/180
	lat2*= pi/180
	long2*= pi/180
    
    Variable dlon = long1 - long2
    Variable dlat = lat1 - lat2
    Variable aa = ((sin(dlat / 2) ^ 2) + cos(lat1) * cos(lat2) * (sin(dlon / 2) ^ 2))
    Variable cc = 2 * atan2(sqrt(aa), sqrt(1 - aa))
    Variable dd = 6371 * cc
    
    // return distance in metres
    return dd * 1000
end

Thirdly there is code to calculate the distance over a ellipsoid using Vincenty's formula, with optional complication of accounting for altitude.

// calculates ellipsoidal distance between two points on earth given coords in decimal degrees
// uses Vincenty's formula + optional approximately correction for altitude above the ellipsoid 
Function EllipsoidalDistanceBetweenTwoPoints(lat1, long1, lat2, long2 [,alt1, alt2])
	Variable lat1, long1, lat2, long2, alt1, alt2
	
	// convert inputs in degrees to radians
	lat1*= pi/180
	long1*= pi/180
	lat2*= pi/180
	long2*= pi/180
	
	// approximately correct for altitude above the ellipsoid 
	variable delta_alt = abs(alt1 - alt2)
	variable avg_alt = (alt1 + alt2) / 2
	
	// correct for errors at exact poles by adjusting 0.6 millimeters
	if (abs(pi/2-abs(lat1)) < 1e-10)
	    lat1 = sign(lat1)*(pi/2-(1e-10))
	endif
	if (abs(pi/2-abs(lat2)) < 1e-10)
	    lat2 = sign(lat2)*(pi/2-(1e-10))
	endif
	
	Variable a = 6378137.0 // equatorial radius in meters 
	Variable f = 1/298.257223563 // ellipsoid flattening 
	Variable b = (1 - f)*a 
	Variable tolerance = 1e-11 // to stop iteration
	
	Variable phi1 = lat1 
	Variable phi2 =lat2
	Variable U1 = atan((1-f)*tan(phi1))
	Variable U2 = atan((1-f)*tan(phi2))
	Variable L1 = long1
	Variable L2 = long2
	Variable L = L2 - L1
	
	Variable lambda_old = L + 0
	Variable t
	Variable sin_sigma
	Variable cos_sigma
	Variable sigma
	Variable sin_alpha
	Variable cos_sq_alpha
	Variable cos_2sigma_m
	Variable C
	Variable lambda_new
	
	do	
		t = (cos(U2)*sin(lambda_old))^2
		t += (cos(U1)*sin(U2) - sin(U1)*cos(U2)*cos(lambda_old))^2
		sin_sigma = t^0.5
		cos_sigma = sin(U1)*sin(U2) + cos(U1)*cos(U2)*cos(lambda_old)
		sigma = atan2(sin_sigma, cos_sigma) 
		
		sin_alpha = cos(U1)*cos(U2)*sin(lambda_old) / sin_sigma
		cos_sq_alpha = 1 - sin_alpha^2
		cos_2sigma_m = cos_sigma - 2*sin(U1)*sin(U2)/cos_sq_alpha
		C = f*cos_sq_alpha*(4 + f*(4-3*cos_sq_alpha))/16
		
		t = sigma + C*sin_sigma*(cos_2sigma_m + C*cos_sigma*(-1 + 2*cos_2sigma_m^2))
		lambda_new = L + (1 - C)*f*sin_alpha*t
		
		 // correct for convergence failure in the case of essentially antipodal points
		if (lambda_new > pi)
			Print "Points are essentially antipodal. Precision may be reduced slightly."
			lambda_new = pi
			break
		endif
			
		if (abs(lambda_new - lambda_old) <= tolerance)
			break
		else
			lambda_old = lambda_new
		endif
	while (1)		
	
	u2 = cos_sq_alpha*((a^2 - b^2)/b^2)
	variable Av = 1 + (u2/16384)*(4096 + u2*(-768+u2*(320 - 175*u2)))
	variable Bv = (u2/1024)*(256 + u2*(-128 + u2*(74 - 47*u2)))
	t = cos_2sigma_m + 0.25*Bv*(cos_sigma*(-1 + 2*cos_2sigma_m^2))
	t -= (Bv/6)*cos_2sigma_m*(-3 + 4*sin_sigma^2)*(-3 + 4*cos_2sigma_m^2)
	Variable delta_sigma = Bv * sin_sigma * t
	Variable s = (b+(avg_alt-(delta_alt/2)))*Av*(sigma - delta_sigma)
	
	// pythagorean theorem
	s = sqrt(s^2 + delta_alt^2)
	
	// in meters	
	return s
end

Some of this is converted from elsewhere and different languages and I take no credit for it.