Warning, /utils/matlab/cs_grid/bk_line/rotate_xy.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit 28aaff14 on 2007-02-05 05:24:33 UTC
28aaff1409 Jean*0001 function [xg1,yg1,xc1,yc1,alp,bet]=rotate_xy(xg0,yg0,xc0,yc0,xa,ya,xb,yb)
0002 %----
0003 % new referential is defined as: equator = great circle AB
0004 % origin of longitude = A , with longitude increasing from A to B.
0005
0006 % find the North pole of the "new" referential.
0007
0008 rad=pi/180;
0009 cla=cos(rad*xa); sla=sin(rad*xa); cpa=cos(rad*ya); spa=sin(rad*ya);
0010 clb=cos(rad*xb); slb=sin(rad*xb); cpb=cos(rad*yb); spb=sin(rad*yb);
0011
0012 x1= cpa*sla*spb-spa*cpb*slb;
0013 x2=-cpa*cla*spb+spa*cpb*clb;
0014 x3=cpa*cpb*sin((xb-xa)*rad);
0015
0016 fac=x1*x1+x2*x2+x3*x3;
0017 bet=asin(x3/sqrt(fac));
0018 alp=atan2(x2,x1);
0019
0020 fprintf('North pole: %8.3f %8.3f\n',alp/rad,bet/rad);
0021
0022 %-----------------
0023
0024 %- rotate grid point coordinate: xg0,yg0 -> xcg,ycg
0025
0026 %- "P" point is on the Eq., on the great circle AB:
0027 % shift longitude, to put North pole at -90.E ("P" point <-> long origin)
0028 xx=xg0-90-alp/rad;
0029
0030 %- rotate North pole, from usual position to new location:
0031 % rotation angle = pi/2 - Beta
0032 cb=cos(pi/2-bet); sb=sin(pi/2-bet);
0033
0034 fprintf('start rotating xG,yG ...');
0035 cx=cos(xx*rad);
0036 sx=sin(xx*rad);
0037 cy=cos(yg0*rad);
0038 sy=sin(yg0*rad);
0039 x1=cy.*cx;
0040 x2=cy.*sx*cb+sy*sb;
0041 x3=-cy.*sx*sb+sy*cb;
0042
0043 yy=asin(x3)/rad;
0044 xx=atan2(x2,x1)/rad;
0045 fprintf(' done\n');
0046
0047 clear x1 x2 x3
0048 cla=cos(rad*xa-pi/2-alp); sla=sin(rad*xa-pi/2-alp); %cpa=cos(rad*ya); spa=sin(rad*ya);
0049 x1=cpa*cla;
0050 x2=cpa*sla*cb+spa*sb;
0051 x3=-cpa*sla*sb+spa*cb;
0052 xa1=atan2(x2,x1)/rad;
0053 ya1=asin(x3)/rad;
0054
0055 clb=cos(rad*xb-pi/2-alp); slb=sin(rad*xb-pi/2-alp); %cpb=cos(rad*yb); spb=sin(rad*yb);
0056 x1=cpb*clb;
0057 x2=cpb*slb*cb+spb*sb;
0058 x3=-cpb*slb*sb+spb*cb;
0059 xb1=atan2(x2,x1)/rad;
0060 yb1=asin(x3)/rad;
0061 fprintf('A pnt, new coord.: %8.3f %8.3f\n',xa1,ya1);
0062 fprintf('B pnt, new coord.: %8.3f %8.3f\n',xb1,yb1);
0063
0064 %- make A point the longitude origin (shift by xa1):
0065 xg1=xx-xa1; yg1=yy;
0066
0067 fprintf('start rotating xC,yC ...');
0068 xx=xc0-90-alp/rad;
0069 cx=cos(xx*rad);
0070 sx=sin(xx*rad);
0071 cy=cos(yc0*rad);
0072 sy=sin(yc0*rad);
0073 x1=cy.*cx;
0074 x2=cy.*sx*cb+sy*sb;
0075 x3=-cy.*sx*sb+sy*cb;
0076
0077 yc1=asin(x3)/rad;
0078 xc1=atan2(x2,x1)/rad - xa1;
0079 fprintf(' done\n');
0080
0081 return
