[1] | 1 | function angle=get_angle(latu,lonu,argu1); |
---|
| 2 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 3 | % |
---|
| 4 | % Compute the grid orientation: angle [radians] |
---|
| 5 | % between XI-axis and the direction to the EAST |
---|
| 6 | % at RHO-points. |
---|
| 7 | % |
---|
| 8 | % lonu longitude of u points |
---|
| 9 | % latu latitude of u points |
---|
| 10 | % argu1: spheroid |
---|
| 11 | % 'clarke66' Clarke 1866 |
---|
| 12 | % 'iau73' IAU 1973 |
---|
| 13 | % 'wgs84' WGS 1984 (default) |
---|
| 14 | % 'sphere' Sphere of radius 6371.0 km |
---|
| 15 | % |
---|
| 16 | % copied from dist.m of the Oceans toolbox |
---|
| 17 | % |
---|
| 18 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 19 | spheroid='wgs84'; |
---|
| 20 | if (nargin >= 3), |
---|
| 21 | spheroid=argu1; |
---|
| 22 | end; |
---|
| 23 | |
---|
| 24 | if (spheroid(1:3)=='sph'), |
---|
| 25 | A = 6371000.0; |
---|
| 26 | B = A; |
---|
| 27 | E = sqrt(A*A-B*B)/A; |
---|
| 28 | EPS= E*E/(1-E*E); |
---|
| 29 | elseif (spheroid(1:3)=='cla'), |
---|
| 30 | A = 6378206.4E0; |
---|
| 31 | B = 6356583.8E0; |
---|
| 32 | E= sqrt(A*A-B*B)/A; |
---|
| 33 | EPS = E*E/(1.-E*E); |
---|
| 34 | elseif(spheroid(1:3)=='iau'), |
---|
| 35 | A = 6378160.e0; |
---|
| 36 | B = 6356774.516E0; |
---|
| 37 | E = sqrt(A*A-B*B)/A; |
---|
| 38 | EPS = E*E/(1.-E*E); |
---|
| 39 | elseif(spheroid(1:3)=='wgs'), |
---|
| 40 | A = 6378137.; |
---|
| 41 | E = 0.081819191; |
---|
| 42 | B = sqrt(A.^2 - (A*E).^2); |
---|
| 43 | EPS= E*E/(1.-E*E); |
---|
| 44 | else |
---|
| 45 | error('dist: Unknown spheroid specified!'); |
---|
| 46 | end; |
---|
| 47 | |
---|
| 48 | latu=latu*pi/180; % convert to radians |
---|
| 49 | lonu=lonu*pi/180; |
---|
| 50 | |
---|
| 51 | latu(latu==0)=eps; % Fixes some nasty 0/0 cases in the |
---|
| 52 | % geodesics stuff |
---|
| 53 | [M,L]=size(latu); |
---|
| 54 | |
---|
| 55 | PHI1=latu(1:M,1:L-1); % endpoints of each segment |
---|
| 56 | XLAM1=lonu(1:M,1:L-1); |
---|
| 57 | PHI2=latu(1:M,2:L); |
---|
| 58 | XLAM2=lonu(1:M,2:L); |
---|
| 59 | |
---|
| 60 | % wiggle lines of constant lat to prevent numerical probs. |
---|
| 61 | PHI2(PHI1==PHI2)=PHI2(PHI1==PHI2)+ 1e-14; |
---|
| 62 | % wiggle lines of constant lon to prevent numerical probs. |
---|
| 63 | XLAM2(XLAM1==XLAM2)=XLAM2(XLAM1==XLAM2)+ 1e-14; |
---|
| 64 | |
---|
| 65 | |
---|
| 66 | % COMPUTE THE RADIUS OF CURVATURE IN THE PRIME VERTICAL FOR |
---|
| 67 | % EACH POINT |
---|
| 68 | |
---|
| 69 | xnu1=A./sqrt(1.0-(E*sin(PHI1)).^2); |
---|
| 70 | xnu2=A./sqrt(1.0-(E*sin(PHI2)).^2); |
---|
| 71 | |
---|
| 72 | % COMPUTE THE AZIMUTHS. azim IS THE AZIMUTH AT POINT 1 |
---|
| 73 | % OF THE NORMAL SECTION CONTAINING THE POINT 2 |
---|
| 74 | |
---|
| 75 | TPSI2=(1.-E*E)*tan(PHI2) + E*E*xnu1.*sin(PHI1)./(xnu2.*cos(PHI2)); |
---|
| 76 | |
---|
| 77 | % SOME FORM OF ANGLE DIFFERENCE COMPUTED HERE?? |
---|
| 78 | |
---|
| 79 | DLAM=XLAM2-XLAM1; |
---|
| 80 | CTA12=(cos(PHI1).*TPSI2 - sin(PHI1).*cos(DLAM))./sin(DLAM); |
---|
| 81 | azim=atan(1./CTA12); |
---|
| 82 | |
---|
| 83 | % GET THE QUADRANT RIGHT |
---|
| 84 | |
---|
| 85 | DLAM2=(abs(DLAM)<pi).*DLAM + (DLAM>=pi).*(-2*pi+DLAM) + ... |
---|
| 86 | (DLAM<=-pi).*(2*pi+DLAM); |
---|
| 87 | azim=azim+(azim<-pi)*2*pi-(azim>=pi)*2*pi; |
---|
| 88 | azim=azim+pi*sign(-azim).*( sign(azim) ~= sign(DLAM2) ); |
---|
| 89 | angle(:,2:L)=(pi/2)-azim(:,1:L-1); |
---|
| 90 | angle(:,1)=angle(:,2); |
---|
| 91 | angle(:,L+1)=angle(:,L); |
---|
| 92 | |
---|
| 93 | return |
---|
| 94 | |
---|