matlab海洋模拟工具箱

function [range,A12,A21]=dist(lat,long,argu1,argu2); % DIST Computes distance and bearing between points on the earth % using various reference spheroids. % % [RANGE,AF,AR]=DIST(LAT,LONG) computes the ranges RANGE between % points specified in the LAT and LONG vectors (decimal degrees with % positive indicating north/east). Forward and reverse bearings % (degrees) are returned in AF, AR. % % [RANGE,GLAT,GLONG]=DIST(LAT,LONG,N) computes N-point geodesics % between successive points. Each successive geodesic occupies % it's own row (N>=2) % % [..]=DIST(...,'ellipsoid') uses the specified ellipsoid % to get distances and bearing. Available ellipsoids are: % % 'clarke66' Clarke 1866 % 'iau73' IAU 1973 % 'wgs84' WGS 1984 % 'sphere' Sphere of radius 6371.0 km % % The default is 'wgs84'. % % Ellipsoid formulas are recommended for distance d<2000 km, % but can be used for longer distances. %Notes: RP (WHOI) 3/Dec/91 % Mostly copied from BDC "dist.f" routine (copied from ....?), but % then wildly modified to bring it in line with Matlab vectorization. % % RP (WHOI) 6/Dec/91 % Feeping Creaturism! - added geodesic computations. This turned % out to be pretty hairy since there were a lot of branch problems % with asin, atan when computing geodesics subtending > 90 degrees % that were ignored in the original code! % RP (WHOI) 15/Jan/91 % Fixed some bothersome special cases, like when computing geodesics % and N=2, or LAT=0... % A Newhall (WHOI) Sep 1997 % modified and fixed a bug found in Matlab version 5 % % NOTE: This routine may interfere with dist that % is supplied with matlab's neural net toolbox. %C GIVEN THE LATITUDES AND LONGITUDES (IN DEG.) IT ASSUMES THE IAU SPHERO %C DEFINED IN THE NOTES ON PAGE 523 OF THE EXPLANATORY SUPPLEMENT TO THE %C AMERICAN EPHEMERIS. %C %C THIS PROGRAM COMPUTES THE DISTANCE ALONG THE NORMAL %C SECTION (IN M.) OF A SPECIFIED REFERENCE SPHEROID GIVEN %C THE GEODETIC LATITUDES AND LONGITUDES OF THE END POINTS %C *** IN DECIMAL DEGREES *** %C %C IT USES ROBBIN'S FORMULA, AS GIVEN BY BOMFORD, GEODESY, %C FOURTH EDITION, P. 122. CORRECT TO ONE PART IN 10**8 %C AT 1600 KM. ERRORS OF 20 M AT 5000 KM. %C %C CHECK: SMITHSONIAN METEOROLOGICAL TABLES, PP. 483 AND 484, %C GIVES LENGTHS OF ONE DEGREE OF LATITUDE AND LONGITUDE %C AS A FUNCTION OF LATITUDE. (SO DOES THE EPHEMERIS ABOVE) %C %C PETER WORCESTER, AS TOLD TO BRUCE CORNUELLE...1983 MAY 27 %C spheroid='wgs84'; geodes=0; if (nargin >= 3), if (isstr(argu1)), spheroid=argu1; else geodes=1; Ngeodes=argu1; if (Ngeodes <2), error('Must have at least 2 points in a goedesic!');end; if (nargin==4), spheroid=argu2; end; end; end; if (spheroid(1:3)=='sph'), A = 6371000.0; B = A; E = sqrt(A*A-B*B)/A; EPS= E*E/(1-E*E); elseif (spheroid(1:3)=='cla'), A = 6378206.4E0; B = 6356583.8E0; E= sqrt(A*A-B*B)/A; EPS = E*E/(1.-E*E); elseif(spheroid(1:3)=='iau'), A = 6378160.e0; B = 6356774.516E0; E = sqrt(A*A-B*B)/A; EPS = E*E/(1.-E*E); elseif(spheroid(1:3)=='wgs'), %c on 9/11/88, Peter Worcester gave me the constants for the %c WGS84 spheroid, and he gave A (semi-major axis), F = (A-B)/A %c (flattening) (where B is the semi-minor axis), and E is the %c eccentricity, E = ( (A**2 - B**2)**.5 )/ A %c the numbers from peter are: A=6378137.; 1/F = 298.257223563 %c E = 0.081819191 A = 6378137.; E = 0.081819191; B = sqrt(A.^2 - (A*E).^2); EPS= E*E/(1.-E*E); else error('dist: Unknown spheroid specified!'); end; NN=max(size(lat)); if (NN ~= max(size(long))), error('dist: Lat, Long vectors of different sizes!'); end if (NN==size(lat)), rowvec=0; % It is easier if things are column vectors, else rowvec=1; end; % but we have to fix things before returning! lat=lat(:)*pi/180; % convert to radians long=long(:)*pi/180; lat(lat==0)=eps*ones(sum(lat==0),1); % Fixes some nasty 0/0 cases in the % geodesics stuff PHI1=lat(1:NN-1); % endpoints of each segment XLAM1=long(1:NN-1); PHI2=lat(2:NN); XLAM2=long(2:NN); % wiggle lines of constant lat to prevent numerical probs. if (any(PHI1==PHI2)), for ii=1:NN-1, if (PHI1(ii)==PHI2(ii)), PHI2(ii)=PHI2(ii)+ 1e-14; end; end; end; % wiggle lines of constant long to prevent numerical probs. if (any(XLAM1==XLAM2)), for ii=1:NN-1, if (XLAM1(ii)==XLAM2(ii)), XLAM2(ii)=XLAM2(ii)+ 1e-14; end; end; end; %C COMPUTE THE RADIUS OF CURVATURE IN THE PRIME VERTICAL FOR %C EACH POINT xnu=A./sqrt(1.0-(E*sin(lat)).^2); xnu1=xnu(1:NN-1); xnu2=xnu(2:NN); %C*** COMPUTE THE AZIMUTHS. A12 (A21) IS THE AZIMUTH AT POINT 1 (2) %C OF THE NORMAL SECTION CONTAINING THE POINT 2 (1) TPSI2=(1.-E*E)*tan(PHI2) + E*E*xnu1.*sin(PHI1)./(xnu2.*cos(PHI2)); PSI2=atan(TPSI2); %C*** SOME FORM OF ANGLE DIFFERENCE COMPUTED HERE?? DPHI2=PHI2-PSI2; DLAM=XLAM2-XLAM1; CTA12=(cos(PHI1).*TPSI2 - sin(PHI1).*cos(DLAM))./sin(DLAM); A12=atan((1.)./CTA12); CTA21P=(sin(PSI2).*cos(DLAM) - cos(PSI2).*tan(PHI1))./sin(DLAM); A21P=atan((1.)./CTA21P); %C GET THE QUADRANT RIGHT DLAM2=(abs(DLAM)<pi).*DLAM + (DLAM>=pi).*(-2*pi+DLAM) + ... (DLAM<=-pi).*(2*pi+DLAM); A12=A12+(A12<-pi)*2*pi-(A12>=pi)*2*pi; A12=A12+pi*sign(-A12).*( sign(A12) ~= sign(DLAM2) ); A21P=A21P+(A21P<-pi)*2*pi-(A21P>=pi)*2*pi; A21P=A21P+pi*sign(-A21P).*( sign(A21P) ~= sign(-DLAM2) ); %%A12*180/pi %%A21P*180/pi SSIG=sin(DLAM).*cos(PSI2)./sin(A12); % At this point we are OK if the angle < 90...but otherwise % we get the wrong branch of asin! % This fudge will correct every case on a sphere, and *almost* % every case on an ellipsoid (wrong hnadling will be when % angle is almost exactly 90 degrees) dd2=[cos(long).*cos(lat) sin(long).*cos(lat) sin(lat)]; dd2=sum((diff(dd2).*diff(dd2))')'; if ( any(abs(dd2-2) < 2*((B-A)/A))^2 ), disp('dist: Warning...point(s) too close to 90 degrees apart'); end; bigbrnch=dd2>2; SIG=asin(SSIG).*(bigbrnch==0) + (pi-asin(SSIG)).*bigbrnch; SSIGC=-sin(DLAM).*cos(PHI1)./sin(A21P); SIGC=asin(SSIGC); A21 = A21P - DPHI2.*sin(A21P).*tan(SIG/2.0); %C COMPUTE RANGE G2=EPS*(sin(PHI1)).^2; G=sqrt(G2); H2=EPS*(cos(PHI1).*cos(A12)).^2; H=sqrt(H2); TERM1=-SIG.*SIG.*H2.*(1.0-H2)/6.0; TERM2=(SIG.^3).*G.*H.*(1.0-2.0*H2)/8.0; TERM3=(SIG.^4).*(H2.*(4.0-7.0*H2)-3.0*G2.*(1.0-7.0*H2))/120.0; TERM4=-(SIG.^5).*G.*H/48.0; range=xnu1.*SIG.*(1.0+TERM1+TERM2+TERM3+TERM4); if (geodes), %c now calculate the locations along the ray path. (for extra accuracy, could %c do it from start to halfway, then from end for the rest, switching from A12 %c to A21... %c started to use Rudoe's formula, page 117 in Bomford...(1980, fourth edition) %c but then went to Clarke's best formula (pg 118) %RP I am doing this twice because this formula doesn't work when we go %past 90 degrees! Ngd1=round(Ngeodes/2); % First time...away from point 1 if (Ngd1>1), wns=ones(1,Ngd1); CP1CA12 = (cos(PHI1).*cos(A12)).^2; R2PRM = -EPS.*CP1CA12; R3PRM = 3.0*EPS.*(1.0-R2PRM).*cos(PHI1).*sin(PHI1).*cos(A12); C1 = R2PRM.*(1.0+R2PRM)/6.0*wns; C2 = R3PRM.*(1.0+3.0*R2PRM)/24.0*wns; R2PRM=R2PRM*wns; R3PRM=R3PRM*wns; %c now have to loop over positions RLRAT = (range./xnu1)*([0:Ngd1-1]/(Ngeodes-1)); THETA = RLRAT.*(1 - (RLRAT.^2).*(C1 - C2.*RLRAT)); C3 = 1.0 - (R2PRM.*(THETA.^2))/2.0 - (R3PRM.*(THETA.^3))/6.0; DSINPSI =(sin(PHI1)*wns).*cos(THETA) + ... ((cos(PHI1).*cos(A12))*wns).*sin(THETA); %try to identify the branch...got to other branch if range> 1/4 circle PSI = asin(DSINPSI); DCOSPSI = cos(PSI); DSINDLA = (sin(A12)*wns).*sin(THETA)