A=pmatequ(T1,T2)

function A=pmatequ(T1,T2)

% function A=pmatequ(T1,T2)
%
% compute matrix of precession  A(i,j) for
% equatorial coordinates from equinox T1 to T2 
% algorithms by Montenbruck/Pfleger: "Astronomy with the Personal Computer"
%
% 05.05.04 M.Penzkofer

rho = pi/180;

A=zeros(3,3);
dT = T2-T1;
ZETA  =  ( (2306.2181+(1.39656-0.000139*T1)*T1)+((0.30188-0.000345*T1)+0.017998*dT)*dT )*dT/3600;
Z     =  ZETA + ( (0.79280+0.000411*T1)+0.000205*dT)*dT*dT/3600;
THETA =  ( (2004.3109-(0.85330+0.000217*T1)*T1)-((0.42665+0.000217*T1)+0.041833*dT)*dT )*dT/3600;
C1=cos(Z*rho);  C2=cos(THETA*rho);  C3=cos(ZETA*rho);
S1=sin(Z*rho);  S2=sin(THETA*rho);  S3=sin(ZETA*rho);
A(1,1)=-S1*S3+C1*C2*C3; A(1,2)=-S1*C3-C1*C2*S3; A(1,3)=-C1*S2;
A(2,1)=+C1*S3+S1*C2*C3; A(2,2)=+C1*C3-S1*C2*S3; A(2,3)=-S1*S2;
A(3,1)=+S2*C3;          A(3,2)=-S2*S3;          A(3,3)=+C2;