j1950toj2000

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   function [ra2000, dec2000] = fk45(r1950, d1950, bepoch)

   function takes data in J1950 coordinates from date bepoch
   (ie. 1970.3) and converts it to J2000.

   r1950,d1950 are in hours and degrees, respectively. as is the output.


   SJCM - copied from some c libraries.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0001 function [ra2000, dec2000] = j1950toj2000(r1950, d1950, bepoch)
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %   function [ra2000, dec2000] = fk45(r1950, d1950, bepoch)
0005 %
0006 %   function takes data in J1950 coordinates from date bepoch
0007 %   (ie. 1970.3) and converts it to J2000.
0008 %
0009 %   r1950,d1950 are in hours and degrees, respectively. as is the output.
0010 %
0011 %
0012 %   SJCM - copied from some c libraries.
0013 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0014 
0015 % definitions
0016 D2PI =  2*pi;
0017 pmf = 100.0 * 60.0 * 60.0 * 360.0 / D2PI;
0018 
0019 a = [ -1.62557e-6,  -0.31919e-6, -0.13843e-6 ];
0020 ad= [1.245e-3,    -1.580e-3,   -0.659e-3 ];
0021 em= [  0.9999256782, -0.0111820611, -0.0048579477 ;
0022        0.0111820610,  0.9999374784, -0.0000271765 ; 
0023        0.0048579479, -0.0000271474,  0.9999881997 ;
0024       -0.000551,     -0.238565,      0.435739     ;
0025        0.238514,     -0.002667,     -0.008541     ;
0026       -0.435623,      0.012254,      0.002117     ];
0027 
0028 % spherical to cartesian
0029 r0(1) = cos(r1950*15*pi/180)*cos(d1950*pi/180);
0030 r0(2) = sin(r1950*15*pi/180)*cos(d1950*pi/180);
0031 r0(3) = sin(d1950*pi/180);
0032 
0033 %  Adjust vector a to give zero proper motion in FK5
0034 w = (bepoch - 1950.0)/pmf;
0035 a1 = a + w*ad;
0036 
0037 % remove e-terms
0038 w = sum(r0.*a1);
0039 v1 = r0 - a1 + w*r0;
0040 %  Convert position vector to Fricke system
0041 for i=1:6
0042   w = 0.0;
0043   for j = 1:3
0044     w = w + em(i,j)*v1(j);
0045   end
0046   v2(i) = w;
0047 end
0048 
0049 % Allow for fictitious proper motion in FK4
0050 mjdDate = 15019.81352 + (bepoch -1900.0)*365.242198781 ;
0051 julianDate = 2000 + (mjdDate -  51544.5 )/365.25;
0052 w = (julianDate - 2000.0)/pmf;
0053 for m=1:3
0054   v2(m) = v2(m) + w*v2(m+3);
0055 end
0056 
0057 %  Revert to spherical coordinates
0058 x = v2(1);
0059 y = v2(2);
0060 z = v2(3);
0061 r = pyth(x,y);
0062 
0063 if(r==0)
0064   w = 0.0;
0065 else
0066   w = atan2(y,x);
0067 end
0068 
0069 if(z==0)
0070   d2000 = 0.0;
0071 else
0072   d2000 = atan2(z,r);
0073 end
0074 
0075 ang = w;
0076 r2000 = mod(w, D2PI);
0077 if(r2000<0)
0078   r2000 = r2000 + D2PI;
0079 end
0080 
0081 % convert to the same input units:
0082 ra2000 = r2000*12/pi;
0083 dec2000 = d2000*180/pi;
0084 
0085 return;