Contents
- Purpose of this page
- Algebraic manipulation of EF = R / (R + h) to solve for R
- Found sample data on NGS web site
- Compute R from manipulated equation
- Compute radii of curvature.
- How sensitive is the determination of EF on the value of R?
- compute difference in ellipsoid for each increment of ten meters
- Summary
% computing radius used in elevation factor determinations % filename: compRfromEF.m % Author: DMMulcare % Date: 16 July 2016 modified 29 July 2016 clear clear all % remove any remaining variables from memory format compact % reduces whitespace when displaying to screen format long g clc % clear screen
Purpose of this page
fprintf('\n\tPrompted by some newgroup discussion about how the radius used in elevation factor (EF) is ') fprintf('\n\tcomputed, I worked through some relevant computations to show ways to compute the radius as well as to ') fprintf('\n\texplore the sensitivity of the EF to changes in R and eh.\n\n')
Prompted by some newgroup discussion about how the radius used in elevation factor (EF) is computed, I worked through some relevant computations to show ways to compute the radius as well as to explore the sensitivity of the EF to changes in R and eh.
Algebraic manipulation of EF = R / (R + h) to solve for R
fprintf('\tFirst we compute a radius (R) given an ellipsoid height(h) and elevation factor (EF)') fprintf('\n\tUsing equation: EF = R / (R + h). \n') fprintf('\n\tRestate equation to solve for R given h and EF. \n') fprintf('\n\t1. Relabel EF to S to avoid confusing it as being two variables S = R / (R + h) \n') fprintf('\t2. Multiply both sides of above by (R + h): S * (R + h) = R \n') fprintf('\t3. Separate terms on left side: (S * R) + (S * h) = R \n') fprintf('\t4. Subtract (S * R) from each side: (S * h) = R - (S * R) \n') fprintf('\t5. Factor right side: (S * h) = R * ( 1 - S) \n') fprintf('\t6. Divide both sides by (1-S): (S * h)/(1 - S) = R \n') fprintf('\t7. Replace S with EF: (EF * h) / (1 - EF) = R \n') fprintf('\n\tLet us verify the manipulation above using real data. \n') %
First we compute a radius (R) given an ellipsoid height(h) and elevation factor (EF) Using equation: EF = R / (R + h). Restate equation to solve for R given h and EF. 1. Relabel EF to S to avoid confusing it as being two variables S = R / (R + h) 2. Multiply both sides of above by (R + h): S * (R + h) = R 3. Separate terms on left side: (S * R) + (S * h) = R 4. Subtract (S * R) from each side: (S * h) = R - (S * R) 5. Factor right side: (S * h) = R * ( 1 - S) 6. Divide both sides by (1-S): (S * h)/(1 - S) = R 7. Replace S with EF: (EF * h) / (1 - EF) = R Let us verify the manipulation above using real data.
Found sample data on NGS web site
extracted data from NGS data sheet for a point atop a lovely vista in Wyoming where I made observations in the early 1990's. Lugging a Trimble 4000SST, two batteries, tripod, cables and more up the last few hundred feet as my two-wheel drive couldn't make it.
fprintf('\n\t MQ0448 DESIGNATION - BLACK BUTTE') fprintf('\n\t MQ0448 STATE/COUNTY- WY/SWEETWATER') fprintf('\n\t MQ0448* NAD 83(2011) POSITION- 41 33 34.71001(N) 108 48 04.50871(W) NO CHECK ') fprintf('\n\t MQ0448* NAD 83(2011) ELLIP HT- 2377.346 (meters) (06/27/12) NO CHECK ') fprintf('\n\t MQ0448; North East Units Scale Factor Converg. ') fprintf('\n\t MQ0448;SPC WYWC - 117,671.878 595,724.815 MT 0.99993772 -0 02 02.4') fprintf('\n\t MQ0448! - Elev Factor x Scale Factor = Combined Factor ') fprintf('\n\t MQ0448!SPC WYWC - 0.99962725 x 0.99993772 = 0.99956500 ')
MQ0448 DESIGNATION - BLACK BUTTE MQ0448 STATE/COUNTY- WY/SWEETWATER MQ0448* NAD 83(2011) POSITION- 41 33 34.71001(N) 108 48 04.50871(W) NO CHECK MQ0448* NAD 83(2011) ELLIP HT- 2377.346 (meters) (06/27/12) NO CHECK MQ0448; North East Units Scale Factor Converg. MQ0448;SPC WYWC - 117,671.878 595,724.815 MT 0.99993772 -0 02 02.4 MQ0448! - Elev Factor x Scale Factor = Combined Factor MQ0448!SPC WYWC - 0.99962725 x 0.99993772 = 0.99956500
Compute R from manipulated equation
fprintf('\n\n\tFirst extract h and EF from the data above.\n') h = 2377.346; EF = 0.99962725; R = (EF * h) / ( 1 - EF); fprintf('\n\tRadius computed from h and EF = %8.2f \n', R); fprintf('\tIs this a reasonable result? \n')
First extract h and EF from the data above. Radius computed from h and EF = 6375479.13 Is this a reasonable result?
Compute radii of curvature.
fprintf('\n\n\tThere are two common radii determinated at a point on the surface of an ellipse.') fprintf('\n\tWe are using GRS80 ellipsoid parameters and our sample point') fprintf('\n\tFirst isthe radius of curvature in the prime vertical (N). \n') fprintf('\tIts equation is: N = a / ( sqrt( 1 - e^2 * sin(latitude)^2 ) \n') fprintf('\n\tNext the radius of curvature in meridian (M). \n') fprintf('\n\tIts equation is: M = (a * (1-e^2)) / (1 - e^2 * sin(lat)^(3/2) ) \n') fprintf('\ta = semi-major axis and e is ellipsoid eccentricity \n') fprintf('\n\n') % Required GRS 80 reference ellipsoid parameters a = 6378137.; % defining parameter finv = 298.257222101; % inverse of flattening eSqr = 6.69438002290341574957494e-3; % this value taken from a web posting by DJMilbert % computed parameters b = a * (1 - (1/finv)); % semi-minor axis e = sqrt(a^2-b^2)/a; % eccentricity ee = e*e; % eccentricity squared diff_eSqr = eSqr - ee; % difference between Matlab and 128-bit value % fprintf('\n\tCompute latitude in radians of the sample point. \n') % tested sample points in WY and FL lat = (41 + 33/60 + 34.71001/3600) * pi/180; % Wyoming %lat = (29 + 5845/10000) * pi/180; % Central Florida N = a / ( sqrt(1 - eSqr * sin(lat)^2) ); M = (a * (1-eSqr)) / (1 - eSqr * sin(lat)^(3/2) ); % Below is equation of radius of ellipse taken from clynchg3c.com Re = (a * sqrt((1-eSqr)^2 * sin(lat)^2 + cos(lat)^2) ) / sqrt(1-eSqr*sin(lat)^2); % Below is equation of average radius taken from Torge text GEODESY printed % page 98 in Google Books image for eBook linked from Wikipedia page Ravg = sqrt(M * N); %This is the Gaussian radius fprintf('\n\tRadius of curvature in the prime vertical at sample point is: %7.1fm \n', N ) fprintf('\t Radius of curvature in the meridian at sample point is: %7.1fm \n', M ) fprintf('\t Radius computed from EF and h is: %7.1fm \n', R ) fprintf('\t Radius of ellipse: %7.1fm \n', Re ) fprintf('\t Average (Gaussian) radius: %7.1fm \n', Ravg) fprintf('\n\tLet us use the solved N and site h to check the site EF value. \n') chkEF = R / (R + h); fprintf('\n\tPublished EF = %2.9f computed EF = %2.9f difference = %1.9f \n', EF, chkEF, EF - chkEF) fprintf('\tComputed R = %7.1fm \tM = %7.1fm \tN = %7.1fm \tRavg = %7.1fm \n',R, M, N, Ravg ) fprintf('\tDifferences R - M = %5.1fm and R - N = %5.1fm and R - Ravg = %7.1fm \n', abs(R-M), abs(R-N), abs(R-Ravg) ) % % shown below are the bounds of the value of N (equator and pole) %Ra = 6378000; % equatorial radius (semi-major axis) %Rpol = 6399593; % N at poles %b = 6356752.31414 % polar radius from b = a(1-f)
There are two common radii determinated at a point on the surface of an ellipse. We are using GRS80 ellipsoid parameters and our sample point First isthe radius of curvature in the prime vertical (N). Its equation is: N = a / ( sqrt( 1 - e^2 * sin(latitude)^2 ) Next the radius of curvature in meridian (M). Its equation is: M = (a * (1-e^2)) / (1 - e^2 * sin(lat)^(3/2) ) a = semi-major axis and e is ellipsoid eccentricity Compute latitude in radians of the sample point. Radius of curvature in the prime vertical at sample point is: 6387553.4m Radius of curvature in the meridian at sample point is: 6358439.1m Radius computed from EF and h is: 6375479.1m Radius of ellipse: 6368769.8m Average (Gaussian) radius: 6372979.6m Let us use the solved N and site h to check the site EF value. Published EF = 0.999627250 computed EF = 0.999627250 difference = 0.000000000 Computed R = 6375479.1m M = 6358439.1m N = 6387553.4m Ravg = 6372979.6m Differences R - M = 17040.0m and R - N = 12074.3m and R - Ravg = 2499.5m
How sensitive is the determination of EF on the value of R?
Use an arbitrary elevation hz = 100;
hz = 100.; Rz = M:1000:N; nz = length(Rz); EFz = zeros(nz,1); fprintf('\n\tRadius in 1KM increments. Right column is elevation factor \n\n') for i=1:nz; EFz(i) = Rz(i)/(Rz(i) + hz); fprintf('\t%8.0f \t\t%+1.8f\n', Rz(i), EFz(i)) end avgDiffEFz = mean(EFz); stdevEFz = std(EFz); fprintf('\n\n\tMean EF is: %+1.8f', avgDiffEFz) fprintf('\n\tStandard deviation: %1.8f\n', stdevEFz) % fprintf('\n\tAnother question arises. How sensitive is the EF to a change in ellipsoid height? \n' ) fprintf('\n\tBelow we determine EF values using our sample site radius R. \n\n') fprintf('\tLeft column is ellipsoid height right column is elevation factor. \n') hx = -40:10:150; n = length(hx); EFx = zeros(nz,1); for i = 1:n; EFx(i) = ( R / (R + hx(i)) ); disp(sprintf('\t%+d \t\t%+1.8f', hx(i), EFx(i) )) end
Radius in 1KM increments. Right column is elevation factor 6358439 +0.99998427 6359439 +0.99998428 6360439 +0.99998428 6361439 +0.99998428 6362439 +0.99998428 6363439 +0.99998429 6364439 +0.99998429 6365439 +0.99998429 6366439 +0.99998429 6367439 +0.99998430 6368439 +0.99998430 6369439 +0.99998430 6370439 +0.99998430 6371439 +0.99998431 6372439 +0.99998431 6373439 +0.99998431 6374439 +0.99998431 6375439 +0.99998432 6376439 +0.99998432 6377439 +0.99998432 6378439 +0.99998432 6379439 +0.99998432 6380439 +0.99998433 6381439 +0.99998433 6382439 +0.99998433 6383439 +0.99998433 6384439 +0.99998434 6385439 +0.99998434 6386439 +0.99998434 6387439 +0.99998434 Mean EF is: +0.99998431 Standard deviation: 0.00000002 Another question arises. How sensitive is the EF to a change in ellipsoid height? Below we determine EF values using our sample site radius R. Left column is ellipsoid height right column is elevation factor. -40 +1.00000627 -30 +1.00000471 -20 +1.00000314 -10 +1.00000157 +0 +1.00000000 +10 +0.99999843 +20 +0.99999686 +30 +0.99999529 +40 +0.99999373 +50 +0.99999216 +60 +0.99999059 +70 +0.99998902 +80 +0.99998745 +90 +0.99998588 +100 +0.99998432 +110 +0.99998275 +120 +0.99998118 +130 +0.99997961 +140 +0.99997804 +150 +0.99997647
compute difference in ellipsoid for each increment of ten meters
diff = zeros(n,1); for i = 1: n-1 diff(i) = EFx(i) - EFx(i+1); % disp(sprintf('\t%+1.9f', diff(i))) % option to display to screen end avgDiffEF = mean(diff); stdevEFdiff = std(diff); fprintf('\n\tAverage difference in EF: %1.5em\n\tstdev: %1.8f \n', avgDiffEF, stdevEFdiff) fprintf('\n\n')
Average difference in EF: 1.49006e-06m stdev: 0.00000035
Summary
fprintf('\n\tFrom the tabulations above we can see that it takes a 4000 meter difference in radius\n\t to change the EF in the eighth decimal place.') fprintf('\n\tWe also see more sensitivity to changes in ellipsoid height where a 10 meter difference \n\t in h leads to a 1.57e-6 change in EF.') fprintf('\n\n')
From the tabulations above we can see that it takes a 4000 meter difference in radius to change the EF in the eighth decimal place. We also see more sensitivity to changes in ellipsoid height where a 10 meter difference in h leads to a 1.57e-6 change in EF.