physsub(PNr,T)

function [MAG,D_EQU,PHI,POSAX,POSSUN] = physsub(PNr,T)

% function [MAG,D_EQU,PHI,POSAX,POSSUN] = physsub(PNr,T)
%
% Physical Ephemeris of the Big Planets and the Sun
% algorithms by Montenbruck/Pfleger: "Astronomy with the Personal Computer"
%
% 05.05.04 M.Penzkofer

% --------------------------------------------------------------------------------------------------------------

% constants
rho = pi/180;
% 1 AU in km
AU      = 149597870.0;
% speed of light (AU/d)
C_LIGHT = 173.14;
% standard epoche J2000
J2000   = 0.0;

% precession matrix (J2000 -> mean equinox of the date) for equatorial coordinates
PMAT = pmatequ(J2000,T);

% equatorial coordinates of the Earth in respect of J2000
[XE,YE,ZE] = position(3,T);
[XE,YE,ZE] = ecl2equ(J2000,XE,YE,ZE);

if (PNr ~= 3)

   	% physical ephemeris of the planets

		% geocentric geometric position of the planet and running time of light (in days)
      [X,Y,Z] = position(PNr,T);
      [X,Y,Z] = ecl2equ(J2000,X,Y,Z);
      DELTA = sqrt ( (X-XE)*(X-XE) + (Y-YE)*(Y-YE) + (Z-ZE)*(Z-ZE) );
      T0 = T - (DELTA/C_LIGHT) / 36525.0;

      % retard position of the planet at the time of light emission
		% (correction for running time of light of 1st order)
      [X,Y,Z] = position(PNr,T0);
      [X,Y,Z] = ecl2equ(J2000,X,Y,Z);

      % geocentric coordinates considering running time of light
      XX = X-XE;
      YY = Y-YE;
      ZZ = Z-ZE;

      % apparent diameter of the equator (in ")
      [R_EQU,F] = shape(PNr);
      D_EQU = 3600.0 * 2.0*asin(R_EQU/(DELTA*AU)) / rho;

      % rectascension and declination of the rotation axis (J2000)
      % and orientation of the zero meridian
      [A0,D0,W1,SENSE] = orient(PNr,1,T0);
      [A0,D0,W2,SENSE] = orient(PNr,2,T0);
      [A0,D0,W3,SENSE] = orient(PNr,3,T0);

      % planetocentric longitude and latitude of the Earth
      [AX,AY,AZ,L1,B,D] = rotation(XX,YY,ZZ,A0,D0,W1,SENSE,F);
      [AX,AY,AZ,L2,B,D] = rotation(XX,YY,ZZ,A0,D0,W2,SENSE,F);
      [AX,AY,AZ,L3,B,D] = rotation(XX,YY,ZZ,A0,D0,W3,SENSE,F);

		% planetocentric longitude and latitude of the Sun
		[AX,AY,AZ,LSUN,BSUN,DSUN] = rotation(X,Y,Z,A0,D0,W1,SENSE,F);

      % illumination parameters and apparent brightness
      [R,DELTA,ELONG,PHI,K] = illum(X,Y,Z,XE,YE,ZE);
      MAG = bright(PNr,R,DELTA,PHI,DSUN,LSUN-L1);

      % position angle of the rotation axis and the Sun
      % in respect of the mean equinox of the date
      PRECART ( PMAT,  X, Y, Z );
      PRECART ( PMAT, XX,YY,ZZ );
      PRECART ( PMAT, AX,AY,AZ );
      POSAX  = POSANG ( XX,YY,ZZ, AX,AY,AZ );
      POSSUN = POSANG ( XX,YY,ZZ, -X,-Y,-Z );

else

		% physical ephemeris of the Sun

		% equatorial heliocentric coordinates of the Earth (J2000)
		% considering the running time of light

		DELTA = sqrt ( XE*XE + YE*YE + ZE*ZE );   
		T0 = T - (DELTA/C_LIGHT) / 36525.0;       
		[XE,YE,ZE] = position(3,T0);
		[XE,YE,ZE] = ecl2equ(J2000,XE,YE,ZE);

		% rectascension and declination of the axis of the Sun (J2000),
		% orientation of the zero meridian and equatorial radius (km)
		A0 = 285.96;  
		D0 =  63.96;
		W1 = 84.11 + 14.1844000 * 36525.0*T;
		W1 = W1/360.0;
		W1 = 360.0*(W1-fix(W1));
		R_EQU = 696000.0;

		% heliographic coordinates of the Earth
		[AX,AY,AZ,L1,B1,B] = rotation(-XE,-YE,-ZE,A0,D0,W1,'RETROGRADE',0.0);

		% position angle of the axis in respect of the mean equinox of the date
		[XE,YE,ZE] = precart(PMAT,XE,YE,ZE);
		[AX,AY,AZ] = precart(PMAT,AX,AY,AZ);
		POSAX = posang(-XE,-YE,-ZE,AX,AY,AZ);

		% position angle of the axis of the Sun in intervall [-180..180] degrees
		if (POSAX > 180.0)
   		POSAX = POSAX-360.0;
		end
      
		% apparent equatorial radius of the Sun (in ")
		D_EQU = 2.0*asin(R_EQU/(DELTA*AU)) / rho;  
      
      % magnitude
      MAG = NaN;
      % phase angle
		PHI = NaN;
      % position angle of the Sun
      POSSUN = NaN;
      
end