% Prabhakar. S
% AE 110 Lab Assignment
% Satellite Tracking Program
% Satellite:  NOAA 12
% Program Name elaz : Elevation-Azimuth being called by the program sts

% Orbital Elements ready for degree to radians conversion!
% Two Line Elements and the Epoch Data are obtained from the GUI Interface

load variables.mat; % Loading the variables from the GUI Interface sts
I2 = yoe;		% Year of epoch
J2 = moe;			% Month of epoch
%h2 = hoe;       % Hour of epoch
%m2 = mioe;       % Minute of epoch
%K2 = doe + hoe/24 + mioe/(24*60);		% Day of epoch
K2=doe;
zero=0; % For printing a single digit month or day preceeded by a zero.
%ut=.08235372;
%frac=.08235372; % fraction of the day from the TLE

% This calculation is to find the Julian Date,
% JD Fortran code written by Fliegel and Van Flandern [1968]
% Handout by Dr. P

JD = (367*I2 - (fix(7*(I2 + fix((J2 + 9)/12))/4)) + fix((275*J2)/9) + K2 + 1721013.5 + frac);

%JD = (K2 - 32075 + round(1461 * (I2 + 4800 + (J2 -14)/12)/4 ...
%+ 367* (J2 - 2 - (J2 - 14)/12*12)/12 - 3*((I2 + 4900 + (J2 - 14)/12)/100)/4)) ;
% fprintf('\n Two Line Elements Data - 28 Dec 2000, Time 11:30:00\n');
% fprintf(' The Julian Date is %f',JD);

% Propagation begins
I = yoe;								% Year
J = moe;								% Month
%h = hoe;   				   		% Hour
%m = mioe;     				 		% Minute
%K = doe + hoe/24 + mioe/(24*60);	% Day
K=doe;

JD0 = (367*I - (fix(7*(I + fix((J + 9)/12))/4)) + fix((275*J)/9) + K + 1721013.5 + frac);

%JD0 = (K - 32075 + round(1461 * (I + 4800 + (J -14)/12)/4 ...
%+ 367 * (J - 2 - (J - 14)/12*12)/12 - 3*((I + 4900 + (J - 14)/12)/100)/4)) ;

I3 = yoe;							% Year
J3 = moe;							% Month
%h3 = hoe;       					% Hour
%m3 = mioe;       				% Minute
%K3 = (doe+pro) + h3/24 + m3/(24*60);	% Day
K3 = (doe+pro);

JD1 = (367*I3 - (fix(7*(I3 + fix((J3 + 9)/12))/4)) + fix((275*J3)/9) + K3 + 1721013.5 + frac);

%JD1 = K3 - 32075 + round(1461 * (I3 + 4800 + (J3 -14)/12)/4 ...
%  + 367 * (J3 - 2 - (J3 - 14)/12*12)/12 - 3*((I3 + 4900 + (J3 - 14)/12)/100)/4);
% Propagation ends here.

% Propagation step Size
step = 5/(24*60); % right now for every 5 minutes
j=1;


fprintf('\n\t\t\t SATELLITE    TRACKING    SYSTEM  \n');
fprintf(' ---------------------------------------------------------------------------');
fprintf('\n				             G r e g o r i a n   D a t e \n');
fprintf(' ---------------------------------------------------------------------------');
fprintf('\n  Julian Date      Azimuth     Elevation  Month/Day/Year Hour:Minute:Seconds\n');
fprintf(' ---------------------------------------------------------------------------');

while JD0<JD1;   
   
% Declaring variables and initializing values 

degtorad = pi/180;		% Radians conversion factor
i = incl*degtorad;		% Inclination
CapOmega = raan*degtorad;		% Right Ascension of the Ascending Node - denoted by symbol Capital Omega
ecc = e;					% Eccentricity
SmallOmega = aop*degtorad;						% Argument of Periapsis - denoted by  small Omega 
Mo = ma * degtorad;		% Initial Mean anomaly
n = revs;	   		% Mean Motion (revs/day)
GP = 3.986e5;				% Gravitational parameter - denoted by nu
a = 24000;					% Semi-major Axis
% a = (GP/n^2)^(1/3);
% fprintf('Semi-major axis is computed to be :\n',a);

Latitude = 0 * degtorad;
Longitude = 280 * degtorad;	
L = ((90*degtorad) - Latitude);			

We = 6.300388097;		% Angular velocity of planet Earth
perturb = 98.9246096453622 * degtorad; % Perturbation due to the Oblateness 
tq = 2448621.5; % JD for 
Qg = perturb + We * (JD0 - tq);
Q = Qg + Longitude;

tp = JD - Mo/n;
M = n * (JD0 - tp);

E = M;				          		% Use M for the first value of E
	f=M-E+ecc*sin(E);			 		% Kepler's Equation
	while abs(f)>0.00000001; 		% Loop for Kepler's equation begins
		f=M-E+ecc*sin(E);			 	%	
		fd=-1+ecc*cos(E);			 	% 
		E=E-f/fd;						%
	end;							 		% loop for Kepler's eqn ends here!


% Computing true anomaly and radius

TA=2*atan((((1+ecc)/(1-ecc))^0.5)*tan(E/2));
    if TA<0;
    	  TA=2*pi+TA;
    end;

r=a*(1-ecc*cos(E)); % radius


% Coordinate transformations

adot = (-(3/2) * n * 1.0827e-3 * (6378.14/a)^2 * cos(i))/(1-ecc^2)^2;
A = CapOmega + adot * (JD0 - JD);
SmallOmega = SmallOmega + adot * (JD0 - JD);

Rseu = [0; 0; 6378.14];
Txyz = [ cos(A) -sin(A) 0; sin(A) cos(A) 0; 0 0 1 ] * [ 1 0 0; 0 cos(i) -sin(i); 0 sin(i) cos(i)] * ...
   [ cos(SmallOmega) -sin(SmallOmega) 0; sin(SmallOmega) cos(SmallOmega) 0; 0 0 1 ];
Tseu = [ cos(L) 0 -sin(L); 0 1 0; sin(L) 0 cos(L)] * [ cos(Q) sin(Q) 0; -sin(Q) cos(Q) 0; 0 0 1];
ruvw = [ (r * cos(TA)); (r * sin(TA)); 0 ];
Pseu = Tseu*(Txyz*ruvw) - Rseu;
P = (Pseu(1)^2 + Pseu(2)^2 + Pseu(3)^2)^0.5;

El = (asin ( Pseu(3)/P ))/degtorad;
Az = (atan2 ( Pseu(2), -Pseu(1)))/degtorad;
	if Az < 0;
   	Az = 360 + Az;
   end;
   
%Loop begins here  %if El>0;
		
array(j,1) = JD0;
array(j,3) = Az;
array(j,4) = El;
%array(j,5) = cat(year,month,day);
j = j+1;
% loop ends here

JD0 = JD0 + step;
			% Gregorian Date from Julian Date 
			% This calculation is valid for any Julian Day Number including negative JDN and 
			% produces a Gregorian date (or possibly a proleptic Gregorian date)

			Z = floor(JD0 - 1721118.5); % JD0 is the Julian Date in propagation 
			R = (JD0 - 1721118.5 - Z);  % R is the fractional part of JD0
			G = (Z - .25);  
			A = (floor(G / 36524.25));  % Calculate the value of A which is the number of full centuries
			B = (A - (A / 4)); % The value of B is this number of days minus a constant
			year = (floor((B+G) / 365.25)); % Calculate the value of Y, the year in a calendar whose years start on March 1
			C = (B + Z - floor(365.25 * year));  % Day count
         month = (fix((5 * C + 456) / 153));  % Month
         UT = (C - fix((153 * month - 457) / 5) + R); % Calculation for UTC
         day = floor(UT);  % Gregorian Day
                    
						if month > 12; 
		      			year = year + 1 ;
		      			month = month - 12; 
                  end; 
                  UT = UT - floor(UT);
         UT = UT*24;
         hr = floor(UT);   %  hour
         UT = UT-floor(UT);
	      UT =UT* 60;
	      min = floor(UT);  % minute
   	   UT =UT- floor(UT);
      	UT =UT* 60;
         secs = round(UT); % seconds 
         if secs==60;
            min=min+1;
            secs=0;
         end;
         % Gregorian Date conversion algorithm ends here
         
         % computations for printintg two digits in case year,month,day,min and secs is < 10
         % (i.e) to print third month, 3 as 03 and day 5 as 05 etc..
         ZERO=int2str(zero);
			MIN=int2str(min); % converting minute into a string
			if min <10;
		  			 MinNew=strcat(ZERO,MIN);
                  min=MinNew;
               else;
                  min=MIN;
               end;
               
         	DAY=int2str(day); % converting DAY into a string
				 if day <10;
		  				 DayNew=strcat(ZERO,DAY);
                   day=DayNew;
                else;
                  day=DAY;
         	 end;
          HR=int2str(hr); % converting Hour into a string
					if hr <10;
		  				 HrNew=strcat(ZERO,HR);
                   hr=HrNew;
               else;
                 	 hr=HR;
		          end;
                
                MONTH=int2str(month); % converting month into a string
						if month <10;
			  				 MonthNew=strcat(ZERO,MONTH);
   	                month=MonthNew;
            	   else
         	         month=MONTH;
                  end;
                  
                 SECS=int2str(secs); % converting seconds into a string
						if secs <10;
			  				 SecsNew=strcat(ZERO,SECS);
   	                secs=SecsNew;
            	   else
         	         secs=SECS;
			         end;

         % computations for printing two digits ends here!
             
         if El>0; 
            fprintf('\n %f | %f | %f |  %s/%s/%4d  |  %s:%s:%s  |',JD0, Az, El, month,day,year, hr,min,secs);
            %plot(El,Az);
         end;   
end;
fprintf('\n ---------------------------------------------------------------------------\n');
fprintf('\n\t\t\t       End of Output\!!! \n\n\n');
