distance between points on the earth
in [Fortran]
|
From: glen herrmannsfeldt on 20 Oct 2009 20:20 analyst41(a)hotmail.com wrote: > Nothing so fancy. Every ZIP code has a latitude and longitude and the > idea is to calculate the distance from one to another. I don't believe that they are known to 16 significant digits. > I am using > FUNCTION ZIPDIST1(RLAT1,RLONG1,RLAT2,RLONG2) > DOUBLE PRECISION RLAT1,RLONG1,RLAT2,RLONG2,ACOSINPUT0,ACOSINPUT A good use for statement functions! DOUBLE PRECISION X SIND(X)=SIN(X/57.2957795130823D0) COSD(X)=COS(X/57.2957795130823D0) > ACOSINPUT0 = & > & SIN(RLAT1/57.2957795130823D0) * SIN(RLAT2/57.2957795130823D0) + COS > (RLAT1/57.2957795130823D0) & > & * COS(RLAT2/57.2957795130823D0) * COS(RLONG2/57.2957795130823D0 - > RLONG1/57.2957795130823D0) > ACOSINPUT = MIN(MAX( -1.0D0,ACOSINPUT0),1.0D0) > ZIPDIST1 = 3958.75586574D0 * acos(acosinput) A few extra significant digits there, too. Is that at sea level? (Remember global warming will change sea level.) > RETURN > END > All the Zipcodes will be in the US. I would say that one percent > accuracy would be more than sufficient - I just want to avoid any > known pitfalls with violent functions such as acos. -- glen
From: Tom Micevski on 23 Oct 2009 23:42 analyst41(a)hotmail.com wrote: > (1) Using Double Precision for latitude and longitude and using double > precision for degrees to radians and radius of the earth (by the way, > I have seen values that differ by about 4 miles) here are some values and references for some common spheroidal models: case ('grs80','gda94') ! Australian ellipsoid a=6378137.00000_rk b=6356752.31414_rk invf=298.257222101_rk f=0.003352810681225_rk case ('wgs84','GPS') ! GPS ellipsoid a=6378137.0_rk b=6356752.3142_rk invf=298.257223563_rk f=one/invf * http://www.icsm.gov.au/icsm/gda/gdatm/gdav2.3.pdf (ICSM, Geocentric Datum of Australia Technical Manual Version 2.3, ISBN 0-9579951-0-5) * http://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf (Appendix A) (NIMA Technical Report TR8350.2, Third Edition, Amendment 1, 3 Jan 2000) > > (2) checking argument of acos to be between -1 and 1 > > Seems to be working OK. > > Are the generic sin,cos and acos functions sufficient? > > Have functions such as DCOS become unnecessary in f95? > > Thanks for any responses. if you're after millimetre accuracy between any 2 points on earth, then you should use the vincenty method: Vincenty (1975) Direct and inverse solutions of geodesics on the ellipsoid with application of nested equations, Survey Review, 23(176), 88-93. some useful links: ! * http://www.ga.gov.au/geodesy/datums/vincenty_inverse.jsp ! Procedure based on: ! * http://www.icsm.gov.au/icsm/gda/gdatm/gdav2.3.pdf (Chapter 4) ! * http://www.movable-type.co.uk/scripts/latlong-vincenty.html ! * http://en.wikipedia.org/wiki/Vincenty's_formulae
From: analyst41 on 24 Oct 2009 06:48 On Oct 23, 11:42 pm, Tom Micevski <n...(a)au-e29b6ec0.invalid> wrote: > analys...(a)hotmail.com wrote: > > (1) Using Double Precision for latitude and longitude and using double > > precision for degrees to radians and radius of the earth (by the way, > > I have seen values that differ by about 4 miles) > > here are some values and references for some common spheroidal models: > > case ('grs80','gda94') ! Australian ellipsoid > a=6378137.00000_rk > b=6356752.31414_rk > invf=298.257222101_rk > f=0.003352810681225_rk > case ('wgs84','GPS') ! GPS ellipsoid > a=6378137.0_rk > b=6356752.3142_rk > invf=298.257223563_rk > f=one/invf > > *http://www.icsm.gov.au/icsm/gda/gdatm/gdav2.3.pdf(ICSM, Geocentric > Datum of Australia Technical Manual Version 2.3, ISBN 0-9579951-0-5) > *http://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf > (Appendix A) (NIMA Technical Report TR8350.2, Third Edition, Amendment > 1, 3 Jan 2000) > > > > > (2) checking argument of acos to be between -1 and 1 > > > Seems to be working OK. > > > Are the generic sin,cos and acos functions sufficient? > > > Have functions such as DCOS become unnecessary in f95? > > > Thanks for any responses. > > if you're after millimetre accuracy between any 2 points on earth, then > you should use the vincenty method: > > Vincenty (1975) Direct and inverse solutions of geodesics on the > ellipsoid with application of nested equations, Survey Review, 23(176), > 88-93. > > some useful links: > ! *http://www.ga.gov.au/geodesy/datums/vincenty_inverse.jsp > ! Procedure based on: > ! *http://www.icsm.gov.au/icsm/gda/gdatm/gdav2.3.pdf(Chapter 4) > ! *http://www.movable-type.co.uk/scripts/latlong-vincenty.html > ! *http://en.wikipedia.org/wiki/Vincenty's_formulae Thanks for the references. Every ZIP code has a latitude and longitude (I guess of some central point in the area covered by the ZIP code). I am using the following code: FUNCTION ZIPDIST1(RLAT1,RLONG1,RLAT2,RLONG2) DOUBLE PRECISION RLAT1,RLONG1,RLAT2,RLONG2,ACOSINPUT0,ACOSINPUT ACOSINPUT0 = & & SIN(RLAT1/57.2957795130823D0) * SIN(RLAT2/57.2957795130823D0) + COS (RLAT1/57.2957795130823D0) & & * COS(RLAT2/57.2957795130823D0) * COS(RLONG2/57.2957795130823D0 - RLONG1/57.2957795130823D0) ACOSINPUT = MIN(MAX( -1.0D0,ACOSINPUT0),1.0D0) ZIPDIST1 = 3958.75586574D0 * acos(acosinput) RETURN END All the Zipcodes will be in the US and 1 percent accuracy will be sufficient. For my purposes, does the Haversine method offer a significant improvement? Thanks.
From: Tom Micevski on 24 Oct 2009 11:33 analyst41(a)hotmail.com wrote: > On Oct 23, 11:42 pm, Tom Micevski <n...(a)au-e29b6ec0.invalid> wrote: >> analys...(a)hotmail.com wrote: >>> (1) Using Double Precision for latitude and longitude and using double >>> precision for degrees to radians and radius of the earth (by the way, >>> I have seen values that differ by about 4 miles) >> here are some values and references for some common spheroidal models: >> >> case ('grs80','gda94') ! Australian ellipsoid >> a=6378137.00000_rk >> b=6356752.31414_rk >> invf=298.257222101_rk >> f=0.003352810681225_rk >> case ('wgs84','GPS') ! GPS ellipsoid >> a=6378137.0_rk >> b=6356752.3142_rk >> invf=298.257223563_rk >> f=one/invf >> >> *http://www.icsm.gov.au/icsm/gda/gdatm/gdav2.3.pdf(ICSM, Geocentric >> Datum of Australia Technical Manual Version 2.3, ISBN 0-9579951-0-5) >> *http://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf >> (Appendix A) (NIMA Technical Report TR8350.2, Third Edition, Amendment >> 1, 3 Jan 2000) >> >> >> >>> (2) checking argument of acos to be between -1 and 1 >>> Seems to be working OK. >>> Are the generic sin,cos and acos functions sufficient? >>> Have functions such as DCOS become unnecessary in f95? >>> Thanks for any responses. >> if you're after millimetre accuracy between any 2 points on earth, then >> you should use the vincenty method: >> >> Vincenty (1975) Direct and inverse solutions of geodesics on the >> ellipsoid with application of nested equations, Survey Review, 23(176), >> 88-93. >> >> some useful links: >> ! *http://www.ga.gov.au/geodesy/datums/vincenty_inverse.jsp >> ! Procedure based on: >> ! *http://www.icsm.gov.au/icsm/gda/gdatm/gdav2.3.pdf(Chapter 4) >> ! *http://www.movable-type.co.uk/scripts/latlong-vincenty.html >> ! *http://en.wikipedia.org/wiki/Vincenty's_formulae > > Thanks for the references. Every ZIP code has a latitude and > longitude (I guess of some central point in the area covered by the > ZIP code). I am using the following code: > > FUNCTION ZIPDIST1(RLAT1,RLONG1,RLAT2,RLONG2) > > > DOUBLE PRECISION RLAT1,RLONG1,RLAT2,RLONG2,ACOSINPUT0,ACOSINPUT > > > ACOSINPUT0 = & > & SIN(RLAT1/57.2957795130823D0) * SIN(RLAT2/57.2957795130823D0) + COS > (RLAT1/57.2957795130823D0) & > & * COS(RLAT2/57.2957795130823D0) * COS(RLONG2/57.2957795130823D0 - > RLONG1/57.2957795130823D0) > > > ACOSINPUT = MIN(MAX( -1.0D0,ACOSINPUT0),1.0D0) > > > ZIPDIST1 = 3958.75586574D0 * acos(acosinput) > > > RETURN > END > > > All the Zipcodes will be in the US and 1 percent accuracy will be > sufficient. For my purposes, does the Haversine method offer a > significant improvement? looks like your using a great circle/spherical trig distance method http://www.ga.gov.au/geodesy/datums/distance.jsp i'm unsure of the accuracy of this method, but the link above gives some indication (although the points are in australia, not usa). it looks like your 1% accuracy will likely be easily met.
First
|
Prev
|
Pages: 1 2 3 4 5 Prev: Integer return values in interfaces Next: Cannot mex with gfortran on Snow Leopard (64bit) |