
; .r rotacia_kupoly_pre_mikroprocesor.pro

; JR, v.1, 2015/02/06

;goto,stat_tmp

; program na vypocet azimutu Slnka pocas jednotlivych dni roka
; pre vypocet parametrov pre rotaciu kupoly

; LSO - Zemepisná poloha: dĺžka - 20° 12' 47.43020", šírka: 49° 11' 42.53495"
lso_long=20.+12./60.+47.43020/3600.
lso_lati=49.+11./60.+42.53495/3600.

; year
yyyy_int=2015
yyyy_str=strcompress(string(yyyy_int),/remove_all)

; calculation of the JD for the beginning of the selected year
;	JDCNV, YR, MN, DAY, HR, JULIAN
; INPUT:
; 	YR = Year, integer scalar or vector
;	MN = Month  integer (1-12) scalar or vector
;	DAY = Day   integer 1-31) scalar or vector 
;	HR  = Hours and fractions of hours of universal time (U.T.), scalar or vector
; OUTPUT:
;	JULIAN = Julian date (double precision) 
;
JDCNV, yyyy_int, 1, 1, 0, jd0

; variables
ndays=366
yr=fltarr(ndays) & mn=yr & day=yr & jd=yr
fod=findgen(24*6)/6./24.
azi=fltarr(ndays,24*6) & hei=azi

; cycle of AZI and HEI calculation for each 10mins in each day

for j=0,ndays-1 do begin ; j - # number of a day in a year DOY

    jd1=jd0+j ; actual JD of this day

    DAYCNV, jd1, yrt, mnt, dayt
    yr(j)=yrt & mn(j)=mnt & day(j)=dayt

    for hh=0,24-1 do begin
        for mm=00,59,10 do begin
            hhmm=hh+mm/60.
            zensun,j,hhmm,lso_lati,lso_long,zenith_angle,azimuth_angle,$
                   solfactor, sunrise=sr,sunset=ss
            ;print,jd1,yrt,mnt,dayt,hhmm,zenith_angle,azimuth_angle
            if j eq 0 then tim(hh*6+mm/10)=hhmm
            azi(j,hh*6+mm/10)=azimuth_angle
            hei(j,hh*6+mm/10)=90.-zenith_angle        
            if azimuth_angle lt 0 then azi(j,hh*6+mm/10)=360.+azimuth_angle
            if solfactor lt 0 then azi(j,hh*6+mm/10)=0.
            if solfactor lt 0 then hei(j,hh*6+mm/10)=0. ;360.+azimuth_angle
        endfor ;mm
    endfor ;hh

endfor ; j -day

stat_tmp:

; summer solstice day = 171
; IDL> p,azi(171,16)-180.     -127.476
; IDL> p,azi(171,112)-180.     127.262
; division of the silstice azimut arc to 255 subarcs of 1 degree width
; arc: no=0 -127 degrees (to E), no=127 0 degrees (S), no 255=128 degrees (to W)

; azi to arc position - positions where the dome orientation is
;                       expected to be adjusted manually 
arc_pos=round(azi-180.)+127
ind=where(hei eq 0.)
arc_pos(ind)=replicate(0.,n_elements(ind))

; calculation of the arc distance made by the Sun in the azimuth
; during a 10min time interval [degrees]
tmp=azi-shift(azi,0,1)
;wtp,(tmp<5.3)>(+1.7)
tmp1=(tmp<5.3)>(+1.7)
indh=where(tmp1 eq 5.3)
indl=where(tmp1 eq 1.7)
tmp2=tmp1
tmp2(indl)=replicate(0.,n_elements(indl))
tmp2(indh)=replicate(0.,n_elements(indh))
arc_dist=shift(tmp2,0,1)        ; arc distance made by the Sun in the azimuth
                                ; during a 10min time interval [degrees]

; calculation of the time interval needed for the Sun to move in
; azimuth for 1 degreee, i.e. for 5.74213
; as 5.74213 cm na kolajnici kupoly zodpoveda 1 uhlovemu stupnu
; (PH,JR,27/1/2015)
stmp=size(arc_dist)                              
time_interval_per_1deg=fltarr(stmp(1),stmp(2))
ind=where(arc_dist ne 0.)
time_interval_per_1deg(ind)=600./arc_dist(ind) ; time interval for Sun motion
                                               ; in azimuth for 1 degree only [s]

; physical [seconds[ to engineering units [0-255] * scale_factor
; IDL> print,max(time_interval_per_1deg)      324.655 s
; IDL> print,min(time_interval_per_1deg)      0.00000 s
scale_factor=max(time_interval_per_1deg)/255. ; 
; IDL> p,scale_factor      1.34375 second/arbitrary 8bit unit
time_interval_per_1deg_eng=round(time_interval_per_1deg/scale_factor)

; output for an microprocessor reading
out=''
for d=0,ndays-1 do begin
    for a=0,24*6-1 do begin
        tmp=time_interval_per_1deg_eng(d,a)  
        tmp1=strcompress(string(tmp),/remove_all)
        if tmp lt 10 then tmp1='00'+tmp1
        if tmp lt 100 and tmp ge 10 then tmp1='0'+tmp1        
        out=out+tmp1
        ;print,d,a,' ',tmp1,' ',out
        ;stop
    endfor
    ;stop
endfor
openw,1,yyyy_str+'_dome_rotation_eng_data.dat'
printf,1,out
close,/all
; format - 1 line of: ndays * (24 hours * 6 10min moments) = 366 * 144
;          = 52704 data when each if of the a3 format, i.e. the file
;          format is 52704a3=a158112 , 158112 bytes
;          example: 000000000000234212....
;          000 means no motion (Sun is below the horizon)
;          234 means 234 * scale_factor [s] after which the dome
;          should move toward W for 1 angular degree, i.e. 5.74213 cm
;
;          scale_factor      1.34375 second/arbitrary 8bit unit
;          file='2015_dome_rotation_eng_data.dat'

print,' End'
print

end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; NAME:
;	JDCNV; PURPOSE:
;	Converts Gregorian dates to Julian days   
;
; EXPLANATION:
;       For IDL versions V5.1 or greater, this procedure is superceded by
;       JULDAY() function in the standard IDL distribution.   Note, however,
;       that prior to V5.1 there wasa bug in JULDAY() that gave answers off
;       by 0.5 days. 
;        
; CALLING SEQUENCE:
;	JDCNV, YR, MN, DAY, HR, JULIAN
;
; INPUTS:
; 	YR = Year, integer scalar or vector
;	MN = Month  integer (1-12) scalar or vector
;	DAY = Day   integer 1-31) scalar or vector 
;	HR  = Hours and fractions of hours of universal time (U.T.), scalar
;              or vector
;		
; OUTPUTS:
;	JULIAN = Julian date (double precision) 
;
; EXAMPLE:
;	To find the Julian Date at 1978 January 1, 0h (U.T.)
;
;	IDL> JDCNV, 1978, 1, 1, 0., JULIAN
;
;	will give JULIAN = 2443509.5
; NOTES:
;	(1) JDCNV will accept vector arguments 
;	(2) JULDATE is an alternate procedure to perform the same function
;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; NAME:
;       SUNPOS
; PURPOSE:
;       To compute the RA and Dec of the Sun at a given date.
;
; CALLING SEQUENCE:
;       SUNPOS, jd, ra, dec, [elong, obliquity, /RADIAN ]
; INPUTS:
;       jd    - The Julian date of the day (and time), scalar or vector
;               usually double precision
; OUTPUTS:
;       ra    - The right ascension of the sun at that date in DEGREES
;               double precision, same number of elements as jd
;       dec   - The declination of the sun at that date in DEGREES
;
; OPTIONAL OUTPUTS:
;       elong - Ecliptic longitude of the sun at that date in DEGREES.
;       obliquity - the obliquity of the ecliptic, in DEGREES
;
; OPTIONAL INPUT KEYWORD:
;       /RADIAN - If this keyword is set and non-zero, then all output variables 
;               are given in Radians rather than Degrees
;
; NOTES:
;       Patrick Wallace (Rutherford Appleton Laboratory, UK) has tested the
;       accuracy of a C adaptation of the sunpos.pro code and found the 
;       following results.   From 1900-2100 SUNPOS  gave 7.3 arcsec maximum 
;       error, 2.6 arcsec RMS.  Over the shorter interval 1950-2050 the figures
;       were 6.4 arcsec max, 2.2 arcsec RMS.  
;
;       The returned RA and Dec are in the given date's equinox.
;
;       Procedure was extensively revised in May 1996, and the new calling
;       sequence is incompatible with the old one.
; METHOD:
;       Uses a truncated version of Newcomb's Sun.    Adapted from the IDL
;       routine SUN_POS by CD Pike, which was adapted from a FORTRAN routine
;       by B. Emerson (RGO).
; EXAMPLE:
;       (1) Find the apparent RA and Dec of the Sun on May 1, 1982
;       
;       IDL> jdcnv, 1982, 5, 1,0 ,jd      ;Find Julian date jd = 2445090.5   
;       IDL> sunpos, jd, ra, dec
;       IDL> print,adstring(ra,dec,2)
;                02 31 32.61  +14 54 34.9
;
;       The Astronomical Almanac gives 02 31 32.58 +14 54 34.9 so the error
;               in SUNPOS for this case is < 0.5".      
;
;       (2) Find the apparent RA and Dec of the Sun for every day in 1997
;
;       IDL> jdcnv, 1997,1,1,0, jd                ;Julian date on Jan 1, 1997
;       IDL> sunpos, jd+ dindgen(365), ra, dec    ;RA and Dec for each day 

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; ROUTINE:      zensun
;
; PURPOSE:      Compute solar position information as a function of
;               geographic coordinates, date and time.
;
; USEAGE:       zensun,day,time,lat,lon,zenith,azimuth,solfac,sunrise,sunset,
;                  local=local
;
; INPUT:
;   day         Julian day (positive scalar or vector)
;               (spring equinox =  80)
;               (summer solstice= 171)
;               (fall equinox   = 266)
;               (winter solstice= 356)
;
;   time        Universal Time in hours (scalar or vector)
;
;   lat         geographic latitude of point on earth's surface (degrees)
;
;   lon         geographic longitude of point on earth's surface (degrees)
;
; OUTPUT:
;
;   zenith      solar zenith angle (degrees)
;
;   azimuth     solar azimuth  (degrees)
;               Azimuth is measured clockwise from due north
;
;   solfac      Solar flux multiplier.  SOLFAC=cosine(ZENITH)/RSUN^2
;               where rsun is the current earth-sun distance in
;               astronomical units.
;
;               NOTE: SOLFAC is negative when the sun is below the horizon
