%%Justin Rice, April 30, 2005
%Cornell University
%jkr25@cornell.edu
%please do not reproduce code with out acknowledgement

clear all


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%user changeable variables
latitude=43*pi/180;
day=60;%day of the year, 0 to 365
numberofpoints=6;%number of data points to produce
startime=10; %starting time of day, 0 to 24
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% end user changeable variables











dt=1;
time=[0:dt:24];
i=1;
while(time(i)<24)
    t = (2 * pi * ((day - 1) / 365.0));
    Declination = (0.322003- 22.971 * cos(t)- 0.357898 * cos(2*t)- 0.14398 * cos(3*t)+ 3.94638 * sin(t)+ 0.019334 * sin(2*t)+ 0.05928 * sin(3*t));
    fDeclination=Declination*pi/180;
    fAltitude = asin(sin(fDeclination) * sin(latitude) + cos(fDeclination) * cos(latitude) * cos((15 * (time(i) - 12)) * (pi / 180)));
    fAzimuth = (acos((cos(latitude) * sin(fDeclination) - cos(fDeclination) * sin(latitude) * cos((15 * (time(i) - 12)) * (pi / 180))) / cos(fAltitude)));
    Altitude(i)=fAltitude*180/pi;
    if(time(i)>12)
        Azimuth(i)=-180+fAzimuth*180/pi;
    else
        Azimuth(i)=180-fAzimuth*180/pi;
    end
    i=i+1;
end
Altitude(i)=Altitude(i-1);
Azimuth(i)=Azimuth(i-1);
result=[time',Altitude',Azimuth'];
%%%%%%%%%%generate hourly data for solarfind
i=startime+1;
j=1;
while(i<=numberofpoints+startime)
    datafortest(j,1)=0;
    datafortest(j,2)=result(i,2)*pi/180;
    datafortest(j,3)=result(i,3)*pi/180;
    i=i+1;
    j=j+1;
end
datafortest
f=0;