// Almanac ~ ~ ~ For a given latitude and longitude and date

// Simon Wheaton-Smith,  Illustrating Shadows Feb 9, 2008
// www.illustratingshadows.com

DegToRad = 2 * 3.1416 / 360;
RadToDeg = 360 / (2 * 3.1416);

// Establish the Julian day  
clc;
mm = input ("Month 1-12:             ");
dd = input ("Day 1-31:               ");
if (mm==1) then; 
  jd = 0+dd; 
end;
if (mm==2) then; 
  jd = 31+dd; 
end;  
if (mm==3) then; 
  jd = 31+28+dd; 
end;  
if (mm==4) then; 
  jd = 31+28+31+dd; 
end;    
if (mm==5) then; 
  jd = 31+28+31+30+dd; 
end;    
if (mm==6) then; 
  jd = 31+28+31+30+31+dd; 
end;  
if (mm==7) then; 
  jd = 31+28+31+30+31+30+dd; 
end;    
if (mm==8) then; 
  jd = 31+28+31+30+31+30+31+dd; 
end;    
if (mm==9) then; 
  jd = 31+28+31+30+31+30+31+31+dd; 
end;    
if (mm==10) then; 
  jd = 31+28+31+30+31+30+31+31+30+dd; 
end;    
if (mm==11) then; 
  jd = 31+28+31+30+31+30+31+31+30+31+dd; 
end;    
if (mm==12) then; 
  jd = 31+28+31+30+31+30+31+31+30+31+30+dd; 
end;    
jdf = jd;

// Get latitude, longitude and reference longitude 
sLat = input ("Latitude [eg 32.75]:    ");
fLat = sLat;
csLat= cos(DegToRad*fLat);

sLng = input ("Longitude [eg 108.2]:   ");
fLng = sLng;

sRef = input ("Ref.Longitude [eg 105]: ");
fRef = sRef;
clng = fLng - fRef;

corm = 4*(fLng-fRef) ; 	
corh = corm / 60 ; 

// Establish the equation of time
eot=-1*(9.84*sin(DegToRad*(2*(360*(jdf-81)/365)))- 7.53*cos(DegToRad*(360*(jdf-81)/365))- 1.5*sin(DegToRad*(360*(jdf-81)/365)))-0.3;

// Establish the sun's declination
decl = RadToDeg*(0.006918 - 0.399912*cos(2*3.1416*(jdf-1)/365) + 0.070257*sin(2*3.1416*(jdf-1)/365)  - 0.006758*cos(2*2*3.1416*(jdf-1)/365) + 0.000907*sin(2*2*3.1416*(jdf-1)/365)  - 0.002697*cos(3*2*3.1416*(jdf-1)/365) + 0.00148*sin(3*2*3.1416*(jdf-1)/365));

// Display useful information
clc;
printf ("*** ALMANAC FOR A DATE AND LOCATION ***\n\n");
printf ("www.illustratingshadows.com\n\n");

printf ("Julian Day:  %3.0f  \n"  , jdf);
printf ("  for Month: %2.0f  \n"  , mm );
printf ("  and Day:   %2.0f  \n\n", dd );

printf ("Latitude:  %5.2f  \n",   fLat);
printf ("Longitude:  %5.2f \n",   fLng);
printf ("Reference:  %5.2f \n\n", fRef);

printf ("EOT mm.m:  %4.2f  ",eot);
if ( eot > 0 ) then;
  printf ("Sun is slow");
end;  
if ( eot < 0 ) then;
  printf ("Sun is fast");
end;

printf ("\n\nDeclination:  %4.1f \n\n", decl);

printf ("Long.Corr mm.m:  %4.2f  \n ", corm);
printf ("Net.Corr mm.m:  %4.2f  \n", (corm+eot));

// Sunrise and set information
//     display rise/set from noon uncorrected

sha = 12-((RadToDeg*(acos( tan(DegToRad*fLat) * tan(DegToRad*decl) )))/15) ;
shh = int(sha);
shm = (sha - int(sha))*60;

printf ("\n[ Sunrise/set & Alt/Azi uncorrected ]\n\n");

// sunrise needs special logic for mins = 0
if (shm == 0) then;
  rizhh = 100*(12-(shh));
else
  rizhh = 100*(12-(shh+1)) + (60-shm);  
end;
 
// sunset needs special logic for mins = 60
sethh = 100*(12+shh)+shm;
if  shm >59 then;
    sethh = 100*(12+shh+1);
end;

printf ( "Rise hhmm:   %3.0f     Set hh.mm   %3.0f     plus ", rizhh, sethh);  
printf ( "net corr: " );
if ( corm+eot >=0 ) then;
 	printf ("+");	          // no need to do if - as implied 	
end; 
printf ( " %4.2f\n\n", (corm+eot) );




for i = 5:12 ;                             
		amh = i;			// am hour 6 to 12	
		alt = RadToDeg*(asin(sin(DegToRad*(decl))*sin(DegToRad*fLat)+cos(DegToRad*decl)*cos(DegToRad*fLat)*cos(DegToRad*(15*(12-(amh))))));
  	 	azi = RadToDeg*(atan(sin(DegToRad*(15*(12-amh)))/(sin(DegToRad*fLat)*cos(DegToRad*(15*(12-amh)))-tan(DegToRad*decl)*cos(DegToRad*fLat))));	  
		printf ("Hour HH: %2i  ", amh);
		printf ("Alt: %5.1f  ", alt);	  		
		printf ("Azi: %5.1f  ", azi);
	    printf ("\n");
end;

printf ("*** ENTER *** to end program\n");
// END


