c shcnox.f c SHC | Navrh ciselniku slunecnich hodin. c SHC | Sundial Design program. c Miroslav Broz (miroslav.broz@email.cz), May 1st 2017 program SLUNECNI_HODINY_CISELNIK c----------------------------------------------------------------------- implicit none c konstanty c constants real*8 pi, pi2, deg, rad parameter(pi = 3.1415926535d0, pi2 = pi/2.d0, deg = pi/180.d0, : rad = 180.d0/pi) c vstupni parametry c input parameters character*80 phistr, lambdastr, dtstr, ddstr, dvzstr real*8 phi, lambda, az, ht, x1, y1, x2, y2, d_nodus, : daz, dht, dah, l_date1 integer gnomon, ndate, t_horiz, t_hodin, t_SEC, t_letni, : t_prodluz, t_delkadne, t_vychod, t_zapad, t_temporal, : t_analema, cast_analemy, t_azimut, t_vyska, t_date1 c interni promenne c internal vars integer i, j, k, m, n, nhour, np, nb, i1st, iflag, nanalema, njd, : nheight, itmp integer iua, iub, iuc, iud, iue, iuf real*8 l, b, dt, dvz, a, h, t, t0, l0, dl, dd, z, : r0(3), a0(3), b0(3), n0(3), phi0(3), r_nodus(3), r0_nodus(3), : r(3), c0(3), Ap, Bp, Cp, An, Bn, Cn, a2(2), r2(2), r_polos2(2), : Ap2(2), Bp2(2), Cp2(2), Dp2(2), hour, minute, second, t_corr, : lambda15, t_vychod2(2), t_zapad2(2), h2(2), eps, h0(3), : rh_nodus(3), tmp, l_arr(25), d, dta, jd, jd0, djd, ra, de, ha, : d2(2), dtt, degree, plocha, az0, ht0 character*80 str,s2 logical lflag, lhour integer n_okraj parameter(n_okraj = 5) real*8 okraj(n_okraj, 2), horiz(n_okraj, 2) integer size parameter(size = 8) character*16 font,empty_file_string,degstr character*5 romannumeral(0:24) c funkce c functions integer length, iargc real*8 reads3, sgn, linterp1, l_delka_dne logical inside_polygon, t_vychod_zapad_chk c data c data data font /'Luxi'/ data degstr /'°'/ data empty_file_string /'("?",/)'/ data romannumeral / : '0', : 'I', : 'II', : 'III', : 'IV', : 'V', : 'VI', : 'VII', : 'VIII', : 'IX', : 'X', : 'XI', : 'XII', : 'XIII', : 'XIV', : 'XV', : 'XVI', : 'XVII', : 'XVIII', : 'XIX', : 'XX', : 'XXI', : 'XXII', : 'XXIII', : 'XXIV' : / c======================================================================= c c nacteni vstupnich parametru ze standardniho vstupu c read input parametres from standard input c if (iargc().eq.0) then i = 0 5 continue i=i+1 read(*,10,err=15,end=15) phistr 10 format(a) if (phistr(1:1).eq.'#') goto 5 i=i+1 read(*,10,err=15,end=15) lambdastr i=i+1 read(*,*,err=15,end=15) az, ht i=i+1 read(*,*,err=15,end=15) x1, y1 i=i+1 read(*,*,err=15,end=15) x2, y2 i=i+1 read(*,*,err=15,end=15) d_nodus i=i+1 read(*,10,err=15,end=15) dtstr i=i+1 read(*,10,err=15,end=15) dvzstr i=i+1 read(*,10,err=15,end=15) ddstr i=i+1 read(*,*,err=15,end=15) daz i=i+1 read(*,*,err=15,end=15) dht i=i+1 read(*,*,err=15,end=15) dah i=i+1 read(*,*,err=15,end=15) ndate i=i+1 read(*,*,err=15,end=15) gnomon i=i+1 read(*,*,err=15,end=15) t_horiz i=i+1 read(*,*,err=15,end=15) t_hodin i=i+1 read(*,*,err=15,end=15) t_SEC i=i+1 read(*,*,err=15,end=15) t_letni i=i+1 read(*,*,err=15,end=15) t_prodluz i=i+1 read(*,*,err=15,end=15) t_delkadne i=i+1 read(*,*,err=15,end=15) t_vychod i=i+1 read(*,*,err=15,end=15) t_zapad i=i+1 read(*,*,err=15,end=15) t_temporal i=i+1 read(*,*,err=15,end=15) t_analema i=i+1 read(*,*,err=15,end=15) cast_analemy i=i+1 read(*,*,err=15,end=15) t_azimut i=i+1 read(*,*,err=15,end=15) t_vyska i=i+1 read(*,*,err=15,end=15) t_date1 i=i+1 read(*,*,err=15,end=15) l_date1 goto 17 15 continue write(*,16) i 16 format('SHC: Chyba pri cteni vstupnich parametru na radku ', : i2,'.') stop 17 continue c c argumenty na prikazovem radku c arguments on the command-line c else if (iargc().eq.30) then i = 1 call getarg(i, phistr) i=i+1 call getarg(i, lambdastr) i=i+1 call getarg(i, str) read(str,*,err=19,end=19) az i=i+1 call getarg(i, str) read(str,*,err=19,end=19) ht i=i+1 call getarg(i, str) read(str,*,err=19,end=19) x1 i=i+1 call getarg(i, str) read(str,*,err=19,end=19) y1 i=i+1 call getarg(i, str) read(str,*,err=19,end=19) x2 i=i+1 call getarg(i, str) read(str,*,err=19,end=19) y2 i=i+1 call getarg(i, str) read(str,*,err=19,end=19) d_nodus i=i+1 call getarg(i, dtstr) i=i+1 call getarg(i, dvzstr) i=i+1 call getarg(i, ddstr) i=i+1 call getarg(i, str) read(str,*,err=19,end=19) daz i=i+1 call getarg(i, str) read(str,*,err=19,end=19) dht i=i+1 call getarg(i, str) read(str,*,err=19,end=19) dah i=i+1 call getarg(i, str) read(str,*,err=19,end=19) ndate i=i+1 call getarg(i, str) read(str,*,err=19,end=19) gnomon i=i+1 call getarg(i, str) read(str,*,err=19,end=19) t_horiz i=i+1 call getarg(i, str) read(str,*,err=19,end=19) t_hodin i=i+1 call getarg(i, str) read(str,*,err=19,end=19) t_SEC i=i+1 call getarg(i, str) read(str,*,err=19,end=19) t_letni i=i+1 call getarg(i, str) read(str,*,err=19,end=19) t_prodluz i=i+1 call getarg(i, str) read(str,*,err=19,end=19) t_delkadne i=i+1 call getarg(i, str) read(str,*,err=19,end=19) t_vychod i=i+1 call getarg(i, str) read(str,*,err=19,end=19) t_zapad i=i+1 call getarg(i, str) read(str,*,err=19,end=19) t_temporal i=i+1 call getarg(i, str) read(str,*,err=19,end=19) t_analema i=i+1 call getarg(i, str) read(str,*,err=19,end=19) cast_analemy i=i+1 call getarg(i, str) read(str,*,err=19,end=19) t_azimut i=i+1 call getarg(i, str) read(str,*,err=19,end=19) t_vyska i=i+1 call getarg(i, str) read(str,*,err=19,end=19) t_date1 i=i+1 call getarg(i, str) read(str,*,err=19,end=19) l_date1 goto 38 19 continue write(*,37) i 37 format('SHC: Chyba pri cteni radkoveho parametru cislo ', : i2,'.') stop 38 continue c vypsani "navodu" c write "help" else write(*,18) 18 format('SHC: Chybny pocet parametru na prikazovem radku. ', : 'Syntaxe parametru:',/,'SHC: phi lambda az ht x1 y1 x2 y2', : ' d_nodus dt dvz dd ndate gnomon t_horiz t_hodin t_SEC', : ' t_letni t_prodluz t_delkadne t_vychod t_zapad t_temporal ', : ' t_analema cast_analemy t_date1 l_date1') stop endif c dekodovani parametru c decode parametres phi = reads3(phistr) lambda = reads3(lambdastr) c konvence takova, aby sedelo cislovani hodin! az = -az*deg ht = ht*deg dt = reads3(dtstr)/pi*180.d0/24.d0 dvz = reads3(dvzstr)/pi*180.d0/24.d0 dd = reads3(ddstr)/pi*180.d0/24.d0 daz = daz*deg dht = dht*deg dah = dah*deg l_date1 = l_date1*deg plocha = abs((x2-x1)*(y2-y1)) c vypis vsechny parametry na stdout pro kontrolu write(*,34) write(*,20) phi*rad write(*,21) lambda*rad write(*,22) -az*rad, ht*rad ! opet minus kvuli zvolene konvenci write(*,23) x1, y1, x2, y2 write(*,33) x2-x1, y2-y1, (x2-x1)/(y2-y1) write(*,24) plocha/1.d4, plocha/100.d0 write(*,25) d_nodus write(*,26) dt, dt*24.d0, dt*1440.d0, dt*86400.d0 write(*,27) dd, dd*24.d0, dd*1440.d0, dd*86400.d0 write(*,28) dvz, dvz*24.d0, dvz*1440.d0, dvz*86400.d0 write(*,32) daz*rad, dht*rad, dah*rad write(*,29) ndate, gnomon, t_horiz, t_hodin, t_SEC, t_letni, : t_prodluz, t_delkadne, t_vychod, t_zapad, t_temporal, : t_analema, cast_analemy, t_azimut, t_vyska, : t_date1, l_date1*rad write(*,255) 34 format('# kontrolni vypis nactenych parametru:') 20 format('# phi = ', f9.5, ' deg') 21 format('# lambda = ', f9.5, ' deg') 22 format('# az = ', f10.5, ' deg ; ht = ', f10.5, ' deg') 23 format('# x1 = ', f9.2, ' cm ; y1 = ', f9.2, ' cm',/, : '# x2 = ', f9.2, ' cm ; y2 = ', f9.2, ' cm') 33 format('# w = ', f9.2, ' cm ; h = ', f9.2, ' cm ;', : ' w/h = ', f8.5) 24 format('# plocha = ', f7.4, ' m^2 = ', f9.2, ' dm^2') 25 format('# d_nodus = ', f9.2, ' cm') 26 format('# dt = ', f8.5, ' day = ', f7.3, ' hour = ', : f8.2, ' min = ', f7.0, ' sec') 27 format('# dd = ', f8.5, ' day = ', f7.3, ' hour = ', : f8.2, ' min = ', f7.0, ' sec') 28 format('# dvz = ', f8.5, ' day = ', f7.3, ' hour = ', : f8.2, ' min = ', f7.0, ' sec') 32 format('# daz = ',f7.3,' deg',/ : '# dht = ',f7.3,' deg',/ : '# dah = ',f7.3,' deg') 29 format( : '# ndate = ',i1,/, : '# gnomon = ',i1,/, : '# t_horiz = ',i1,/, : '# t_hodin = ',i1,/, : '# t_SEC = ',i1,/, : '# t_letni = ',i1,/, : '# t_prodluz = ',i1,/, : '# t_delkadne = ',i1,/, : '# t_vychod = ',i1,/, : '# t_zapad = ',i1,/, : '# t_temporal = ',i1,/, : '# t_analema = ',i1,/, : '# cast_analemy = ',i1,/, : '# t_azimut = ',i1,/, : '# t_vyska = ',i1,/, : '# t_date1 = ',i1,/, : '# l_date1 = ',f7.3,' deg') 255 format('#') c======================================================================= c vystupni unity c input units iua = 10 iub = 20 iuc = 30 iud = 40 iue = 50 iuf = 60 c zaokrouhleni lambda na 15 deg, pokud ne SEC! lambda15 = int(lambda*rad/15.d0+0.5)*15.d0*deg c vystupni soubor s casovou korekci open(unit=iua,file='korekce.dat',status='unknown') write(iua,35) 35 format('# casova korekce na zemepisnou delku: minuty & sekundy') if (t_SEC.ne.1) then t_corr = (lambda15-lambda)/pi*12.d0 call hhms(abs(t_corr)+0.5d0/3600.d0,hour,minute,second) write(iua,36) int(sgn(t_corr)*minute), int(second) 36 format(i4,1x,i2) lambda = lambda15 else write(iua,36) 0, 0 endif close(iua) c casova korekce vzhledem k UT pro popisky hodinovych rysek t_corr = lambda15*rad/15.d0 if (t_letni.eq.1) then t_corr = t_corr + 1.d0 endif c plocha ciselniku r = r0 + A*a0 + B*b0 c 2 vektory v plose a0(1) = 0 a0(2) = -1 a0(3) = 0 b0(1) = -1 b0(2) = 0 b0(3) = 0 c normala n0(1) = 0 n0(2) = 0 n0(3) = 1 c otoceni ciselniku (dle azimutu a vysky normaly) z = ht-pi2 call rotat0(a0,az,z) call rotat0(b0,az,z) call rotat0(n0,az,z) write(*,251) 251 format('# rovnice roviny ciselniku: r = r0 + a0*x + b0*y + n0*z') write(*,250) 'a0',a0(1), a0(2), a0(3) 250 format('# ',a,' = [ ',f14.9,' , ',f14.9,' , ',f14.9,']') write(*,250) 'b0',b0(1), b0(2), b0(3) write(*,250) 'n0',n0(1), n0(2), n0(3) write(*,255) c smer rotacni osy Zeme phi0(1) = -cos(phi) phi0(2) = 0.d0 phi0(3) = sin(phi) c datovy soubor pro kresleni ukazatele open(unit=iua,file='nodus.dat',status='unknown') write(iua,45) 45 format('# souradnice paty ukazatele a nodu (x, y, z) [cm]', : ' vzhledem k rovine ciselniku') c zde budou ukladany vsechny popisky grafu open(unit=iue,file='popisky.plt',status='unknown') write(iue,46) 46 format('# popisky vyznacnych bodu pro Gnuplot (PostScript)') c a zde se vypisi vsechny textove popisky open(unit=iuf,file='popisky.dat',status='unknown') write(iuf,49) 49 format('# vsechny popisky rysek ciselniku (viz tez popisky.ps)') c======================================================================= c c ukazatel POLOS c POLOS indicator c if (gnomon.eq.0) then c pocatek, pata ciselniku r0(1) = 0.d0 r0(2) = 0.d0 r0(3) = 0.d0 c totez v rovine ciselniku (odtud vychazeji hodinove cary!) r_polos2(1) = 0.d0 r_polos2(2) = 0.d0 c ukazatel smeruje vzdy NAD plochu (skalarni soucin vektoru phi . n) tmp = 0.d0 do i = 1, 3 tmp = tmp + phi0(i)*n0(i) enddo if (tmp.lt.0) then d_nodus = -d_nodus endif c poloha nodu do i = 1, 3 r_nodus(i) = r0(i) + d_nodus*phi0(i) enddo c prumet nodu do roviny ciselniku call intersect_plane_line(r0, a0, b0, r_nodus, n0, An, Bn, Cn, : lflag) do i = 1, 3 r0_nodus(i) = r_nodus(i) + Cn*n0(i) enddo c kolmy prumetu nodu (souradnice nodu v rovine ciselniku a jeho vyska nad plochou) write(iua,48) 0.d0, 0.d0, 0.d0 48 format(f9.2,1x,f9.2,1x,f9.2) write(iua,48) An, Bn, -Cn write(iue,47) 'P', 0.d0,0.d0 write(iue,47) 'N', An,Bn write(iue,47) 'Ns', -Cn,Bn 47 format('set label " ',a,'" at ',g16.8,',',g16.8,' left front') c======================================================================= c c ukazatel GNOMON c GNOMON indicator c else c pocatek, pata gnomonu r0_nodus(1) = 0.d0 r0_nodus(2) = 0.d0 r0_nodus(3) = 0.d0 c poloha nodu do i = 1, 3 r_nodus(i) = r0_nodus(i) + d_nodus*n0(i) enddo c nalezeni paty polosu jako pruseciku roviny ciselniku a r_nodus + Cn*phi call intersect_plane_line(r0_nodus, a0, b0, r_nodus, phi0, : An, Bn, Cn, lflag) do i = 1, 3 r0(i) = r0_nodus(i) + An*a0(i) + Bn*b0(i) enddo c prepis r0 opet na pocatek! r0(1) = 0.d0 r0(2) = 0.d0 r0(3) = 0.d0 c potrebuji pouze souradnice paty v rovine ciselniku (odkud vychazeji hodinove rysky) r_polos2(1) = An r_polos2(2) = Bn c nakresleni gnomonu v Gnuplotu write(iua,48) 0.d0, 0.d0, 0.d0 write(iua,48) 0.d0, 0.d0, d_nodus write(iue,47) 'N', 0.d0,0.d0 write(iue,47) 'Ns', d_nodus,0.d0 write(*,252) 252 format('# souradnice vyznacnych bodu pro vypocet kulisoveho ', : 'ukazatele:') write(*,253) 'P', An, Bn, 0.d0 253 format('# ',a,' = [ ',g16.8,' , ',g16.8,' , ',g16.8,']') write(*,253) 'N', 0.d0, 0.d0, d_nodus write(*,253) 'Ns', 0.d0, 0.d0, 0.d0 write(*,255) endif close(iua) c======================================================================= c c OBDELNIK CISELNIKU jako datovy soubor pro gnuplot c RECTANGLE OF THE DIAL as a gnuplot data file c open(unit=iua,file='ciselnik.dat',status='unknown') write(iua,41) 41 format('# souradnice bodu okraje ciselniku: x [cm] & y [cm]') write(iua,40) x1,y1 write(iua,40) x1,y2 write(iua,40) x2,y2 write(iua,40) x2,y1 write(iua,40) x1,y1 40 format(f8.2,1x,f8.2) close(iua) open(unit=iua,file='ciselnik.plt',status='unknown') write(iua,51) 51 format('# nastaveni rozsahu os x a y v Gnuplotu', : ' = rozmer ciselniku [cm]') write(iua,50) x1,x2,y1,y2 50 format('x1 = ', f8.2,/,'x2 = ',f8.2,/ : 'y1 = ', f8.2,/,'y2 = ',f8.2,/, : 'set xr [x1:x2]',/,'set yr [y1:y2]') close(iua) c definice polygonu okraj(1,1) = x1 okraj(1,2) = y1 okraj(2,1) = x2 okraj(2,2) = y1 okraj(3,1) = x2 okraj(3,2) = y2 okraj(4,1) = x1 okraj(4,2) = y2 okraj(5,1) = x1 okraj(5,2) = y1 c======================================================================= c c OBRAZ HORIZONTU c PROJECTION OF THE HORIZON c c vodorovny prumet nodu do roviny ciselniku h0(1) = 1.d0 h0(2) = 0.d0 h0(3) = 0.d0 call intersect_plane_line(r0, a0, b0, r_nodus, h0, An, Bn, Cn, : lflag) if (lflag) then do i = 1, 3 rh_nodus(i) = r_nodus(i) + Cn*h0(i) enddo c pruseciky vodorovne cary v ciselniku s okrajem r2(1) = An r2(2) = Bn a2(1) = 1.d0 a2(2) = 0.d0 np = 2 call intersect_polygon_line_2D(okraj,n_okraj,r2,a2,2,Ap2,Bp2,np) else An = 0.d0 Bn = 0.d0 Cn = 0.d0 np = 0 endif open(unit=iua,file='horiz.dat',status='unknown') open(unit=iub,file='horiz_rysky.dat',status='unknown') write(iua,140) write(iub,141) 140 format('# vodorovny prumet nodu do roviny ciselniku', : ' a vzdálenost od nodu: x [cm] & y [cm] & z [cm]') 141 format('# pruseciky horizontu s okrajem ciselniku: ', : 'x [cm] & y [cm]') write(iua,*) An, Bn, Cn do k = 1, np write(iub,*) Ap2(k), Bp2(k) enddo write(iub,empty_file_string) close(iua) close(iub) do i = 1, n_okraj do j = 1, 2 horiz(i,j) = okraj(i,j) enddo enddo c modifikace okraje ciselniku (pro hodiny od V/Z slunce) if ((np.eq.2).and.(t_horiz.eq.0)) then eps = 1.d-3 c okraj(3,2) = Bp2(1)+eps c okraj(4,2) = Bp2(2)+eps horiz(3,2) = Bp2(1)+eps horiz(4,2) = Bp2(2)+eps endif c======================================================================= c c HODINOVE ZNACKY A RYSKY (PRO ZIMNI A LETNI SLUNOVRAT) c HOUR LABELS AND MARKS (FOR THE SUMMER AND WINTER SOLSTICE) c l0 = 90.d0*deg dl = 180.d0*deg b = 0.d0*deg nhour = int(1.d0/dt+0.5d0) eps = 1.d-4 if (t_hodin.eq.0) then nhour = 0 endif open(unit=iua,file='hodin.dat',status='unknown') open(unit=iub,file='hodin_prumet.dat',status='unknown') open(unit=iuc,file='hodin_rysky.dat',status='unknown') write(iua,60) write(iub,61) write(iuc,62) 60 format('# polohy stredniho slunce: hodina UT [h] & azimut [deg]', : ' & vyska [deg]') 61 format('# hodinove rysky: hodina [h] & azimut [deg] & ', : 'vyska [deg] stredniho slunce',/, : '# & x [cm] & y [cm] poloha prumetu nodu na ciselniku', : ' & vzdalenost prumetu',/,'# od nodu [cm]') 62 format('# pruseciky hodinovych rysek s okrajem ciselniku: ', : 'hodina [h] & x [cm] & y [cm]') c polohy stredniho slunce pro vsechny hodiny (pocitany zvlast) do j = 1, 2 l = l0 + (j-1)*dl write(iua,30) l*rad write(iub,30) l*rad do i = 0, nhour t = i*dt*24.d0 call stredni_slunce_ah(l,b,t,lambda,phi,a,h) write(iua,64) t, a*rad, h*rad 64 format(f7.4,1x,f8.3,1x,f7.3) enddo enddo c prumety znacek do roviny ciselniku do i = 0, nhour t = i*dt*24.d0 c kontrola, zda ALESPON NEKDY behem roku slunce sviti nad ciselnikem c check, if AT LEAST SOMETIMES during the year, the Sun is above the dial lhour = .false. n = 180 do j = 0, n l = l0 + j*180.0d0*deg/n call stredni_slunce_ah(l,b,t,lambda,phi,a,h) c0(1) = cos(a)*cos(h) c0(2) = sin(a)*cos(h) c0(3) = sin(h) call intersect_plane_line(r0,a0,b0,r_nodus,c0, : Ap2(1),Bp2(1),Cp2(1),lflag) c velke kladne Cp2 je zde kvuli rysce pro 6 h na J hodinach c large positive Cp2 is here due to the line for 6th hour on a southern dial if (((Cp2(1).lt.0d0).or.(Cp2(1).gt.1.d8)).and. : ((t_horiz.eq.1).or.(h.gt.-eps))) lhour = .true. enddo c znacky pro letni/zimni slunovrat c labels for summer/winter solstice do j = 1, 2 l = l0 + (j-1)*dl call stredni_slunce_ah(l,b,t,lambda,phi,a,h) eps = 1.d-4 h = h + eps c0(1) = cos(a)*cos(h) c0(2) = sin(a)*cos(h) c0(3) = sin(h) h = h - eps call intersect_plane_line(r0,a0,b0,r_nodus,c0, : Ap2(j),Bp2(j),Cp2(j),lflag) a2(j) = a h2(j) = h enddo c pruseciky hodinovych rysek s okrajem ciselniku c intersections of hour lines with dial border if (lhour) then do k = 1, 2 if ((Cp2(k).lt.0.d0).and. : ((t_horiz.eq.1).or.(h2(k).gt.-eps))) then write(iub,65) t + t_corr,a2(k)*rad,h2(k)*rad, : Ap2(k),Bp2(k),Cp2(k) 65 format(f7.4,1x,f8.3,1x,f7.3,1x,g12.4,1x,g12.4,1x,g12.4) endif enddo write(iub,*) call hodinove_rysky(Ap2,Bp2,Cp2,okraj,horiz,n_okraj, : t_prodluz,r_polos2,nb) do k = 1, nb Ap = Ap2(k) Bp = Bp2(k) c popisky jen na ciselniku a MIMO patu polosu (pri protazeni k) c labels only inside the dial and OUTSIDE the polos base (in case of extension to) if ((inside_polygon(okraj,n_okraj,Ap,Bp,1)).and. : (sqrt((Ap-r_polos2(1))**2+(Bp-r_polos2(2))**2).gt.eps)) : then c ciselne popisky hodin pro Gnuplot call hhms24(t + t_corr + 0.5d0/3600.d0,hour,minute,second) if (int(minute).eq.0) then write(iue,90) int(hour),Ap,Bp,font(1:length(font)),size 90 format('set label "',i2,' " at ',g16.8,',',g16.8, : ' right rotate font "',a,',',i3,'" front') write(iuf,92) int(hour) 92 format(i2) else write(iue,91) int(hour),int(minute),Ap,Bp, : font(1:length(font)),size 91 format('set label "',i2,'^{',i2,'} " at ', : g16.8,',',g16.8,' right rotate font "',a,',',i3, : '" front') write(iuf,93) int(hour),int(minute) 93 format(i2,':',i2) endif endif ! inside_polygon write(iuc,70) t + t_corr,Ap,Bp 70 format(f7.4,1x,g12.4,1x,g12.4) enddo write(iuc,*) endif enddo c a missing data character due to gnuplot 4.0.0 write(iua,empty_file_string) write(iub,empty_file_string) write(iuc,empty_file_string) close(iua) close(iub) close(iuc) close(iud) write(iuf,*) c======================================================================= c c DATOVE KRIVKY c DATE CURVES c c klasicke krivky po 30 (nebo 90) stupnich ekliptikalni delky if (t_delkadne.eq.0) then if (ndate.eq.0) then dl = 0.d0 else if (ndate.eq.3) then dl = 90.d0 else if (ndate.eq.7) then dl = 30.d0 else write(*,*) 'SHC: Nepodporovany pocet datovych krivek', : '(0, 3 nebo 7).' stop endif l0 = 90.d0*deg dl = dl*deg do i = 1, ndate l_arr(i) = l0 + (i-1)*dl enddo c krivky pro celociselne delky dne else l_arr(1) = 90.d0*deg l_arr(2) = 270.d0*deg j = 2 do i = 1, 24 call delka_dne_l(i*1.d0, phi, l, iflag) if (iflag.eq.0) then j = j + 1 l_arr(j) = l + 180.d0*deg endif enddo ndate = j endif if (t_date1.eq.1) then ndate = ndate + 1 l_arr(ndate) = l_date1 endif nhour = int(1.d0/dd+0.5d0) open(unit=iua,file='datum.dat',status='unknown') open(unit=iub,file='datum_prumet.dat',status='unknown') write(iua,80) write(iub,81) 80 format('# polohy stredniho slunce: cas [h] & azimut [deg] & ', : 'vyska [deg] & cislo datove cary') 81 format('# datove krivky: hodina [h] & azimut [deg] & ', : 'vyska [deg] stredniho slunce',/, : '# & x [cm] & y [cm] poloha prumetu nodu na ciselniku', : ' & vzdalenost prumetu',/,'# od nodu [cm]') do j = 1, ndate l = l_arr(j) d = l_delka_dne(phi, l, b) i1st = 0 write(iua,30) l*rad write(iub,30) l*rad 30 format('# l = ', f7.3, ' deg') do i = 0, nhour t = i*dd*24.d0 c pro lambda = 0! call stredni_slunce_ah(l,b,t,0.d0,phi,a,h) write(iua,31) t, a*rad, h*rad, j 31 format(f7.4,1x,f8.3,1x,f7.3,1x,i2) c prumety datovych car do roviny ciselniku c0(1) = cos(a)*cos(h) c0(2) = sin(a)*cos(h) c0(3) = sin(h) call intersect_plane_line(r0,a0,b0,r_nodus,c0,Ap,Bp,Cp,lflag) c do k = 1, 3 c r(k) = r_nodus(k) + Cp*c0(k) c enddo if ((Cp.lt.0.d0).and.((h.gt.-eps).or.(t_horiz.eq.1))) then write(iub,65) t,a*rad,h*rad,Ap,Bp,Cp c popisek datove krivky v Gnuplotu if ((Ap.gt.x1).and.(Ap.lt.x2).and.(Bp.gt.y1).and.(Bp.lt.y2) : .and.(i1st.eq.0)) then i1st = 1 if (t_delkadne.eq.0) then str = ' ' else call hhms(d + 0.5d0/3600.d0,hour,minute,second) if (int(minute).eq.0) then write(str,97) int(hour) 97 format(i2) else write(str,98) int(hour),int(minute) 98 format(i2,':',i2) endif call str_trim(str,' ') s2 = '(' // str(1:length(str)) // ')' str = s2 endif write(iue,95) int(l*rad+0.5d0/3600.d0), : degstr(1:length(degstr)),str(1:length(str)), : Ap,Bp,font(1:length(font)),size 95 format('set label " ',i3,a,' ',a,'" at ', : g16.8,',',g16.8,' left font "',a,',',i3,'" front') write(iuf,96) int(l*rad),char(176) 96 format(i3,a1) if (t_delkadne.ne.0) then write(iuf,10) str(1:length(str)) endif endif endif enddo write(iua,*) write(iub,*) enddo write(iua,empty_file_string) write(iub,empty_file_string) close(iua) close(iub) write(iuf,*) c======================================================================= c c RYSKY HODIN POCITANYCH OD VYCHODU SLUNCE c HOURS CALCULATED FROM THE SUNRISE c open(unit=iub,file='vychod_prumet.dat',status='unknown') open(unit=iuc,file='vychod_rysky.dat',status='unknown') write(iub,121) write(iuc,122) 121 format('# rysky od V slunce: hodina UT [h] & azimut [deg] & ', : 'vyska [deg] stredniho slunce',/, : '# & x [cm] & y [cm] poloha prumetu nodu na ciselniku', : ' & vzdalenost prumetu',/,'# od nodu [cm]') 122 format('# pruseciky rysek od V slunce s okrajem ciselniku: ', : 'hodina UT [h] & x [cm] & y [cm]') c casy vychodu/zapadu slunce pro letni/zimni slunovrat c times of sunrise/sunset for summer/winter solstice l0 = 90.d0*deg dl = 180.d0*deg b = 0.d0*deg nhour = int(1.d0/dvz+0.5d0) eps = 1.d-3 do j = 1, 2 l = l0 + (j-1)*dl write(iub,30) l*rad call t_vychod_zapad(l,b,lambda,phi,t_vychod2(j),t_zapad2(j)) call hhms(t_vychod2(j) + 0.5d0/60.d0,hour,minute,second) write(iub,110) int(hour),int(minute) 110 format('# t_vychod = ',i2,' h ',i2,' min') enddo if (t_vychod.eq.0) then n = -1 else n = nhour endif do i = 0, n-1 do j = 1, 2 t0 = i*dvz*24.d0 t = t0 + t_vychod2(j) l = l0 + (j-1)*dl c zde uz musi but lambda = 0, protoze jsem delku zapocital jiz do t c lambda = 0 has to be here, because the longitude is already included in t call stredni_slunce_ah(l,b,t,0.d0,phi,a,h) c prumety znacek do roviny ciselniku c0(1) = cos(a)*cos(h) c0(2) = sin(a)*cos(h) c0(3) = sin(h) call intersect_plane_line(r0,a0,b0,r_nodus,c0, : Ap2(j),Bp2(j),Cp2(j),lflag) a2(j) = a h2(j) = h enddo c kontrola, zda alespon NEKDY behem roku slunce sviti nad ciselnikem c check, if at least SOMETIMES during the year, the Sun is above the dial lhour = t_vychod_zapad_chk(t0,lambda,phi,r0,a0,b0,r_nodus, : t_horiz,1) if (lhour) then do k = 1, 2 write(iub,65) t0,a2(k)*rad,h2(k)*rad,Ap2(k),Bp2(k),Cp2(k) enddo write(iub,*) c spolecny kod pro V/Z hodiny call hodinove_rysky(Ap2,Bp2,Cp2,okraj,horiz,n_okraj, : 0,r_polos2,nb) do k = 1, nb write(iuc,70) t0,Ap2(k),Bp2(k) enddo write(iuc,*) c ciselne popisky hodin pro Gnuplot call hhms(t0 + 0.5d0/3600.d0,hour,minute,second) c popisek nekde na care tmp = 0.25d0 Ap = linterp1(Ap2(1),Ap2(2),tmp,0) Bp = linterp1(Bp2(1),Bp2(2),tmp,0) c popisky jen na ciselniku if ((nb.gt.0).and.(inside_polygon(okraj,n_okraj,Ap,Bp,1))) : then str = romannumeral(int(hour)) if (int(minute).eq.0) then write(iue,130) str(1:length(str)),Ap,Bp, : font(1:length(font)),size 130 format('set label "',a,' " at ',g16.8,',',g16.8, : ' center rotate font "',a,',',i3,'" front') write(iuf,132) str(1:length(str)) 132 format(a) else write(iue,131) str(1:length(str)),int(minute),Ap,Bp, : font(1:length(font)),size 131 format('set label "',a,'^{',i2,'} " at ', : g16.8,',',g16.8,' center rotate font "',a,',',i3, : '" front') write(iuf,133) str(1:length(str)),int(minute) 133 format(a,':',i2) endif endif endif enddo write(iub,empty_file_string) write(iuc,empty_file_string) close(iub) close(iuc) write(iuf,*) c======================================================================= c c RYSKY HODIN POCITANYCH OD ZAPADU SLUNCE PREDCHOZIHO DNE c HOURS CALCULATED FROM THE SUNSET OF THE PREVIOUS DAY c open(unit=iub,file='zapad_prumet.dat',status='unknown') open(unit=iuc,file='zapad_rysky.dat',status='unknown') write(iub,151) write(iuc,152) 151 format('# rysky od Z slunce: hodina UT [h] & azimut [deg] & ', : 'vyska [deg] stredniho slunce',/, : '# & x [cm] & y [cm] poloha prumetu nodu na ciselniku', : ' & vzdalenost prumetu',/,'# od nodu [cm]') 152 format('# pruseciky rysek od Z slunce s okrajem ciselniku: ', : 'hodina UT [h] & x [cm] & y [cm]') do j = 1, 2 l = l0 + (j-1)*dl write(iub,30) l*rad call hhms(t_zapad2(j) + 0.5d0/60.d0,hour,minute,second) write(iub,115) int(hour),int(minute) 115 format('# t_zapad = ',i2,' h ',i2,' min') enddo if (t_zapad.eq.0) then n = -1 else n = nhour endif do i = 1, n do j = 1, 2 t0 = i*dvz*24.d0 t = t0 + t_zapad2(j) l = l0 + (j-1)*dl c zde uz musi but lambda = 0, protoze jsem delku zapocital jiz do t c lambda = 0 has to be here, because the longitude is already included in t call stredni_slunce_ah(l,b,t,0.d0,phi,a,h) c prumety znacek do roviny ciselniku c projections of points into the dial plane c0(1) = cos(a)*cos(h) c0(2) = sin(a)*cos(h) c0(3) = sin(h) call intersect_plane_line(r0,a0,b0,r_nodus,c0, : Ap2(j),Bp2(j),Cp2(j),lflag) a2(j) = a h2(j) = h enddo c kontrola, zda ALESPON NEKDY behem roku slunce sviti nad ciselnikem c check, if AT LEAST SOMETIMES during the year, the Sun is above the dial lhour = t_vychod_zapad_chk(t0,lambda,phi,r0,a0,b0,r_nodus, : t_horiz,0) if (lhour) then do k = 1, 2 write(iub,65) t0,a2(k)*rad,h2(k)*rad,Ap2(k),Bp2(k),Cp2(k) enddo write(iub,*) call hodinove_rysky(Ap2,Bp2,Cp2,okraj,horiz,n_okraj, : 0,r_polos2,nb) do k = 1, nb write(iuc,70) t0,Ap2(k),Bp2(k) enddo write(iuc,*) c ciselne popisky hodin pro Gnuplot call hhms(t0 + 0.5d0/3600.d0,hour,minute,second) c popisek na care tmp = 0.75d0 Ap = linterp1(Ap2(1),Ap2(2),tmp,0) Bp = linterp1(Bp2(1),Bp2(2),tmp,0) c popisky jen na ciselniku if ((nb.gt.0).and.(inside_polygon(okraj,n_okraj,Ap,Bp,1))) : then str = romannumeral(int(hour)) if (int(minute).eq.0) then write(iue,160) str(1:length(str)),Ap,Bp, : font(1:length(font)),size 160 format('set label "',a,' " at ',g16.8,',',g16.8, : ' center rotate font "',a,',',i3,'" front') write(iuf,132) str(1:length(str)) else write(iue,161) str(1:length(str)),int(minute),Ap,Bp, : font(1:length(font)),size 161 format('set label "',a,'^{',i2,'} " at ', : g16.8,',',g16.8,' center rotate font "',a,',',i3, : '" front') write(iuf,133) str(1:length(str)),int(minute) endif endif endif enddo write(iub,empty_file_string) write(iuc,empty_file_string) close(iub) close(iuc) write(iuf,*) c======================================================================= c c ANALEMA c ANALEMMA c dta = dt if (t_analema.eq.0) then nanalema = 0 else if (t_analema.eq.1) then nanalema = 1 t0 = 12.d0 - t_corr else nanalema = int(1.d0/dta+0.5d0) t0 = 0.d0 endif c referencni rok je 2000 (resp. zimni slunovrat 1999) jd0 = 2451544.5d0 - 11.d0 djd = 1.d0 njd = int(366.d0/djd+0.5d0) if (cast_analemy.eq.1) then njd = njd/2 jd0 = jd0 + njd else if (cast_analemy.eq.2) then njd = njd/2 endif open(unit=iua,file='analema.dat',status='unknown') open(unit=iub,file='analema_prumet.dat',status='unknown') write(iua,170) write(iub,171) 170 format('# polohy praveho slunce: julianske datum [JD] & ', : 'azimut [deg] & vyska [deg] & cislo analemy') 171 format('# analemy: julianske datum [JD] & azimut [deg] & ', : 'vyska [deg] stredniho slunce',/, : '# & x [cm] & y [cm] poloha prumetu nodu na ciselniku', : ' & vzdalenost prumetu',/,'# od nodu [cm]') do j = 1, nanalema t = t0 + (j-1)*dta*24.d0 i1st = 0 write(iua,180) t write(iub,180) t 180 format('# t = ', f7.3, ' h') do i = 0, njd jd = jd0 + i*djd + t/24.d0 call sunah(jd,lambda,phi,1,0,0,a,h,ra,de,ha,l,b) write(iua,181) jd, a*rad, h*rad, j 181 format(f13.5,1x,f8.3,1x,f7.3,1x,i4) c prumety analem do roviny ciselniku c0(1) = cos(a)*cos(h) c0(2) = sin(a)*cos(h) c0(3) = sin(h) call intersect_plane_line(r0,a0,b0,r_nodus,c0,Ap,Bp,Cp,lflag) if ((Cp.lt.0.d0).and.((h.gt.-eps).or.(t_horiz.eq.1))) then write(iub,182) jd,a*rad,h*rad,Ap,Bp,Cp 182 format(f13.5,1x,f8.3,1x,f7.3,1x,g12.4,1x,g12.4,1x,g12.4) endif enddo write(iua,*) write(iub,*) enddo write(iua,empty_file_string) write(iub,empty_file_string) close(iua) close(iub) c======================================================================= c c TEMPORALNI HODINY c TEMPORAL HOURS c nhour = 12 dtt = dvz l0 = 90.d0*deg dl = 180.d0*deg b = 0.d0*deg open(unit=iub,file='temporal_prumet.dat',status='unknown') open(unit=iuc,file='temporal_rysky.dat',status='unknown') write(iub,191) write(iuc,192) 191 format('# rysky temporalnich hodin: hodina UT [h] & ', : 'azimut [deg] & vyska [deg] stredniho slunce',/, : '# & x [cm] & y [cm] poloha prumetu nodu na ciselniku', : ' & vzdalenost prumetu',/,'# od nodu [cm]') 192 format('# pruseciky rysek temporalnich hodin s okrajem', : ' ciselniku: hodina UT [h] & x [cm] & y [cm]') do j = 1, 2 l = l0 + (j-1)*dl write(iub,30) l*rad d2(j) = l_delka_dne(phi,l,b) call hhms(d2(j) + 0.5d0/60.d0,hour,minute,second) write(iub,195) int(hour),int(minute) 195 format('# delka_dne = ',i2,' h ',i2,' min') d2(j) = d2(j)/nhour*(dtt*24.d0) enddo if (t_temporal.eq.0) then n = -1 else n = int(nhour/(dtt*24.d0)+0.5d0) endif do i = 0, n do j = 1, 2 t0 = dtt*i*24.d0 + 6.d0 if (t0.gt.12.d0+eps) then t0 = t0 - 12.d0 endif t = t_vychod2(j) + i*d2(j) l = l0 + (j-1)*dl call stredni_slunce_ah(l,b,t,0.d0,phi,a,h) c prumety znacek do roviny ciselniku c0(1) = cos(a)*cos(h) c0(2) = sin(a)*cos(h) c0(3) = sin(h) call intersect_plane_line(r0,a0,b0,r_nodus,c0, : Ap2(j),Bp2(j),Cp2(j),lflag) a2(j) = a h2(j) = h enddo c alespon jedna znacka nad obzorem a ciselnikem if (((Cp2(1).lt.0.d0).and.(h2(1).gt.-eps)).or. : ((Cp2(2).lt.0.d0).and.(h2(2).gt.-eps))) then do k = 1, 2 write(iub,65) t0,a2(k)*rad,h2(k)*rad,Ap2(k),Bp2(k),Cp2(k) enddo write(iub,*) call hodinove_rysky(Ap2,Bp2,Cp2,okraj,horiz,n_okraj, : 0,r_polos2,nb) do k = 1, nb write(iuc,70) t0,Ap2(k),Bp2(k) enddo write(iuc,*) c ciselne popisky hodin pro Gnuplot call hhms(t0 + 0.5d0/3600.d0,hour,minute,second) c popisek na care tmp = 0.5d0 Ap = linterp1(Ap2(1),Ap2(2),tmp,iflag) Bp = linterp1(Bp2(1),Bp2(2),tmp,iflag) c popisky jen na ciselniku if ((nb.gt.0).and.(inside_polygon(okraj,n_okraj,Ap,Bp,1))) : then if (int(minute).eq.0) then write(iue,200) int(hour),Ap,Bp,font(1:length(font)),size 200 format('set label "',i2,' " at ',g16.8,',',g16.8, : ' center font "',a,',',i3,'" front') write(iuf,92) int(hour) else write(iue,201) int(hour),int(minute),Ap,Bp, : font(1:length(font)),size 201 format('set label "',i2,'^{',i2,'} " at ', : g16.8,',',g16.8,' center font "',a,',',i3,'" front') write(iuf,93) int(hour),int(minute) endif endif endif enddo write(iub,empty_file_string) write(iuc,empty_file_string) close(iub) close(iuc) write(iuf,*) c======================================================================= c c AZIMUTALNI USECKY c AZIMUTH ABSCISSAS c open(unit=iub,file='azimut_prumet.dat',status='unknown') open(unit=iuc,file='azimut_rysky.dat',status='unknown') write(iub,210) write(iuc,211) 210 format('# rysky azimutalnich usecek: ', : 'azimut [deg] & vyska [deg] & x [cm] & y [cm] ',/, : '# poloha prumetu nodu na ciselniku & ', : 'vzdalenost prumetu od nodu [cm]') 211 format('# pruseciky azimutalnich usecek s okrajem', : ' ciselniku: azimut [deg] & x [cm] & y [cm]') if (t_azimut.eq.0) then n = -1 else n = int(360.d0*deg/daz+0.5d0) endif eps = 1.d-2 do i = 0, n a = 0.d0 + i*daz c moznost zobrazeni neskutecnych usecek nad urovni horizontu c optionally, display unreal lines above the horizon do m = 0, t_horiz do j = 1, 2 c h nabyva postupne hodnot 0, 90, 0, -90 stupnu c h should sequentially have values 0, 90, 0, and -90 degrees h = (j-1)*90.d0*deg * (1-2*m) c prumety znacek do roviny ciselniku c projections of marks into the dial plane c posun vysky tak, aby nebyla celociselna (nulova) c shift the height slightly from integer (zero) values h = h - eps * ((j-1)*2-1) c0(1) = cos(a)*cos(h) c0(2) = sin(a)*cos(h) c0(3) = sin(h) c h = h + eps * ((j-1)*2-1) call intersect_plane_line(r0,a0,b0,r_nodus,c0, : Ap2(j),Bp2(j),Cp2(j),lflag) a2(j) = a h2(j) = h enddo c alespon jedna znacka nad obzorem a ciselnikem c at least one mark above the horizon and the dial if (((Cp2(1).lt.0.d0).and.(h2(1).gt.-eps)).or. : ((Cp2(2).lt.0.d0).and.(h2(2).gt.-eps))) then do k = 1, 2 write(iub,220) a2(k)*rad,h2(k)*rad,Ap2(k),Bp2(k),Cp2(k) 220 format(f8.3,1x,f7.3,1x,g12.4,1x,g12.4,1x,g12.4) enddo write(iub,*) call hodinove_rysky(Ap2,Bp2,Cp2,okraj,horiz,n_okraj, : 0,r_polos2,nb) do k = 1, nb write(iuc,221) a*rad,Ap2(k),Bp2(k) 221 format(f9.4,1x,g12.4,1x,g12.4) enddo write(iuc,*) c ciselne popisky hodin pro Gnuplot c popisky POUZE pro realne azimutalni usecky (pro slunce nad obzorem) if (m.eq.0) then c popisek vzdy na OKRAJI ciselniku if (inside_polygon(okraj,n_okraj,Ap2(1),Bp2(1),-1)) then tmp = 1.d0 else tmp = 0.d0 endif Ap = linterp1(Ap2(1),Ap2(2),tmp,iflag) Bp = linterp1(Bp2(1),Bp2(2),tmp,iflag) c popisky jen na ciselniku a jen pro cary pod horizontem if ((nb.gt.0).and.inside_polygon(okraj,n_okraj,Ap,Bp,1)) : then degree = a*rad if (degree.gt.180.d0+eps) degree = degree - 360.d0 write(str,222) int(degree+sgn(degree)*0.5d0),degstr 222 format(i4,a) write(iue,223) str(1:length(str)),Ap,Bp, : font(1:length(font)),size 223 format('set label "',a,' " at ',g16.8,',',g16.8, : ' right rotate font "',a,',',i3,'" front') call str_replace(str,'.',',') write(iuf,10) str(1:length(str)) endif endif endif ! nad/pod horizontem/ciselnikem enddo ! t_horiz enddo write(iub,empty_file_string) write(iuc,empty_file_string) close(iub) close(iuc) write(iuf,*) c======================================================================= c c VYSKOVE KRIVKY c ALTITUDE CURVES c open(unit=iub,file='vyska_prumet.dat',status='unknown') write(iub,230) 230 format('# vyskove krivky: ', : 'azimut [deg] & vyska [deg] & x [cm] & y [cm] ',/, : '# poloha prumetu nodu na ciselniku & ', : 'vzdalenost prumetu od nodu [cm]') if (t_vyska.eq.0) then n = -1 else if (t_horiz.eq.0) then n = int(90.d0*deg/dht+0.5d0) ht0 = 0.d0*deg else n = int(180.d0*deg/dht+0.5d0) ht0 = -90.d0*deg endif endif nheight = int(360.d0*deg/dah+0.5d0) do i = 0, n h = ht0 + i*dht i1st = 0 write(iub, 231) h*rad 231 format('# h = ',f7.3,' deg') if (abs(abs(h)-90.d0*deg).gt.eps) then m = nheight az0 = az + 180.d0*deg else m = 0 az0 = 0.d0 endif do j = 0, m a = az0 - j*dah c prumety znacek do roviny ciselniku c projections of marks into the dial plane c0(1) = cos(a)*cos(h) c0(2) = sin(a)*cos(h) c0(3) = sin(h) call intersect_plane_line(r0,a0,b0,r_nodus,c0, : Ap,Bp,Cp,lflag) if (Cp.lt.0.d0) then write(iub,220) -a*rad,h*rad,Ap,Bp,Cp c popisek vyskove krivky v Gnuplotu c label of the height curve in Gnuplot if ((Ap.gt.x1).and.(Ap.lt.x2).and.(Bp.gt.y1).and.(Bp.lt.y2) : .and.(i1st.eq.0)) then i1st = 1 write(str,222) int(h*rad+sgn(h)*0.5d0),degstr write(iue,224) str(1:length(str)),Ap,Bp, : font(1:length(font)),size 224 format('set label "',a,' " at ',g16.8,',',g16.8, : ' left font "',a,',',i3,'" front') call str_replace(str,'.',',') write(iuf,10) str(1:length(str)) endif endif enddo write(iub,*) enddo write(iub,empty_file_string) close(iub) c write(iuf,*) c Nakonec zavri soubory s popisky... c Finally close files with labels... close(iue) close(iuf) stop end c======================================================================= c----------------------------------------------------------------------- c c Pruseciky hodinovych (i jinych) rysek s okrajem ciselniku. c Intersection of hour (and other) lines with dial border. c c Ap2, Bp2, Cp2 ... souradnice prumetu x, y, z pro 2 rysky c okraj ... polygon popisujici okraj ciselniku c horiz ... okraj oriznuty horizontem c t_prodluz ... rysky jako 0 usecky, .ge.2 protazeni k okraji, c .eq.1.or.3 k pate polosu c r_polos2 ... souradnice paty polosu v rovine ciselniku c nb ... vysledny pocet pruseciku subroutine hodinove_rysky(Ap2,Bp2,Cp2,okraj,horiz,n_okraj, : t_prodluz,r_polos2,nb) implicit none integer n_okraj, t_prodluz, nb, i1st real*8 Ap2(2), Bp2(2), Cp2(2), okraj(n_okraj), horiz(n_okraj), : r_polos2(2) integer np real*8 C2(2), D2(2), r2(2), a2(2), tmp c functions logical inside_polygon nb = 0 ! na zacatku zadne hodinove rysky nemame... c c nejprve pripad jednoduche USECKY (slunce sviti NAD ciselnikem v zime i v lete) c if ((Cp2(1).lt.0).and.(Cp2(2).lt.0)) then c je alespon jeden z bodu uvnitr ciselniku? if ((inside_polygon(horiz,n_okraj,Ap2(1),Bp2(1),-1)) : .or.(inside_polygon(horiz,n_okraj,Ap2(2),Bp2(2),-1))) then i1st = 1 else i1st = 0 endif c jako pocatecni bod ten, ktery je uvnitr ciselniku! if ((i1st.eq.1).and. : (.not.(inside_polygon(horiz,n_okraj,Ap2(1),Bp2(1),-1)))) then tmp = Ap2(1) Ap2(1) = Ap2(2) Ap2(2) = tmp tmp = Bp2(1) Bp2(1) = Bp2(2) Bp2(2) = tmp endif r2(1) = Ap2(1) r2(2) = Bp2(1) a2(1) = Ap2(2) - r2(1) a2(2) = Bp2(2) - r2(2) np = 2 call intersect_polygon_line_2D(horiz,n_okraj,r2,a2,0, : C2,D2,np) if (np.eq.2) then Ap2(1) = C2(1) Bp2(1) = D2(1) Ap2(2) = C2(2) Bp2(2) = D2(2) nb = 2 elseif (np.eq.1) then Ap2(1) = r2(1) Bp2(1) = r2(2) Ap2(2) = C2(1) Bp2(2) = D2(1) nb = 2 elseif (i1st.eq.1) then ! zadne pruseciky, ale oba body jsou uvnitr nb = 2 endif c jako prvni bod ten blize paty polosu if (((Ap2(1)-r_polos2(1))**2+(Bp2(1)-r_polos2(2))**2).gt. : ((Ap2(2)-r_polos2(1))**2+(Bp2(2)-r_polos2(2))**2)) then tmp = Ap2(1) Ap2(1) = Ap2(2) Ap2(2) = tmp tmp = Bp2(1) Bp2(1) = Bp2(2) Bp2(2) = tmp tmp = Cp2(1) Cp2(1) = Cp2(2) Cp2(2) = tmp endif c c prodlouzit rysky k okraji ciselniku? c if (t_prodluz.ge.2) then r2(1) = Ap2(1) r2(2) = Bp2(1) a2(1) = Ap2(2) - r2(1) a2(2) = Bp2(2) - r2(2) np = 2 call intersect_polygon_line_2D(okraj,n_okraj,r2,a2,1, : C2,D2,np) if (np.eq.2) then Ap2(1) = C2(1) Bp2(1) = D2(1) Ap2(2) = C2(2) Bp2(2) = D2(2) nb = 2 elseif (np.eq.1) then Ap2(1) = r2(1) Bp2(1) = r2(2) Ap2(2) = C2(1) Bp2(2) = D2(1) nb = 2 endif endif ! t_prodluz c c prodlouzit rysky k pate polosu? c if ((t_prodluz.eq.1).or.(t_prodluz.eq.3)) then r2(1) = r_polos2(1) r2(2) = r_polos2(2) a2(1) = Ap2(2) - r2(1) a2(2) = Bp2(2) - r2(2) np = 2 call intersect_polygon_line_2D(okraj,n_okraj,r2,a2,1, : C2,D2,np) if (np.eq.2) then Ap2(1) = C2(1) Bp2(1) = D2(1) Ap2(2) = C2(2) Bp2(2) = D2(2) nb = 2 elseif (np.eq.1) then Ap2(1) = r2(1) Bp2(1) = r2(2) Ap2(2) = C2(1) Bp2(2) = D2(1) nb = 2 endif endif ! t_prodluz c c POLOPRIMKY - pripad, kdyz je slunce v zime nebo v lete jiz POD ciselnikem c else c jako prvni bod ten, pro nejz je slunce nad ciselnikem if (Cp2(1).gt.0.d0) then tmp = Ap2(1) Ap2(1) = Ap2(2) Ap2(2) = tmp tmp = Bp2(1) Bp2(1) = Bp2(2) Bp2(2) = tmp tmp = Cp2(1) Cp2(1) = Cp2(2) Cp2(2) = tmp endif r2(1) = Ap2(1) r2(2) = Bp2(1) a2(1) = r2(1) - Ap2(2) a2(2) = r2(2) - Bp2(2) np = 2 call intersect_polygon_line_2D(horiz,n_okraj,r2,a2,1, : C2,D2,np) if (np.eq.2) then Ap2(1) = C2(1) Bp2(1) = D2(1) Ap2(2) = C2(2) Bp2(2) = D2(2) nb = 2 elseif (np.eq.1) then Ap2(1) = r2(1) Bp2(1) = r2(2) Ap2(2) = C2(1) Bp2(2) = D2(1) nb = 2 endif c k okraji ciselniku neni treba nic prodluzovat! :-) c zbyva tedy prodlouzeni rysek k pate polosu if ((t_prodluz.eq.1).or.(t_prodluz.eq.3)) then r2(1) = r_polos2(1) r2(2) = r_polos2(2) a2(1) = Ap2(1) - r2(1) a2(2) = Bp2(1) - r2(2) np = 2 call intersect_polygon_line_2D(okraj,n_okraj,r2,a2,1, : C2,D2,np) if (np.eq.2) then Ap2(1) = C2(1) Bp2(1) = D2(1) Ap2(2) = C2(2) Bp2(2) = D2(2) nb = 2 elseif (np.eq.1) then Ap2(1) = r2(1) Bp2(1) = r2(2) Ap2(2) = C2(1) Bp2(2) = D2(1) nb = 2 endif endif ! t_prodluz endif ! Cp2 return end c----------------------------------------------------------------------- c c Prusecik(y) usecky (urcene bodem a vektorem: rp = r + a) c s polygonem (2D pole souradnic) v TEZE rovine. Pro iflag.eq.1 hleda c pruseciky uvnitr poloprimky; iflag.eq.2 naopak protina primku. c subroutine intersect_polygon_line_2D(body,n,r,a,iflag,Ap,Bp,np) implicit none integer n,iflag,np real*8 body(n,2),r(2),a(2),Ap(np),Bp(np) integer i,j real*8 r0(2),b(2),Cp,Dp,x,y logical prusecik c constants real*8 eps,infty parameter(infty = 1.d38) c functions real*8 sgn,plusminus eps = 1.d-38 do i = 1, np Ap(i) = r(1) Bp(i) = r(2) enddo j = 0 Cp = -1.d0 Dp = -1.d0 do i = 1, n-1 r0(1) = body(i,1) r0(2) = body(i,2) b(1) = body(i+1,1) - r0(1) b(2) = body(i+1,2) - r0(2) x = (r(1)-r0(1))*a(2) - (r(2)-r0(2))*a(1) y = b(1)*a(2) - b(2)*a(1) if (abs(y).gt.eps) then Cp = x/y else Cp = plusminus(x)*plusminus(y)*infty endif x = r0(1)-r(1)+Cp*b(1) y = a(1) if (abs(y).gt.eps) then Dp = x/y else Dp = plusminus(x)*plusminus(y)*infty endif if (j.eq.np) then np = j return endif eps = 1.d-4 ! tato hodnota nesmi byt prilis velika prusecik = .false. c primka if (iflag.eq.2) then if ((Cp.ge.0.d0).and.(Cp.le.1.d0)) then prusecik = .true. endif c poloprimka elseif (iflag.eq.1) then if ((Cp.ge.0.d0).and.(Cp.le.1.d0).and.(Dp.gt.(-eps))) then prusecik = .true. endif c usecka else if ((Cp.ge.0.d0).and.(Cp.le.1.d0).and. : (Dp.gt.(-eps)).and.(Dp.lt.(1.d0+eps))) then prusecik = .true. endif endif if (prusecik) then j = j + 1 Ap(j) = r0(1) + Cp*b(1) Bp(j) = r0(2) + Cp*b(2) endif enddo np = j return end c----------------------------------------------------------------------- c c Test, zda je bod uvnitr polygonu (zatim pouze obdelniku v zakladni c poloze!). sig [+-1] urcuje, zda hranice patri "dovnitr" ci "ven". c logical function inside_polygon(body, n, x, y, sig) implicit none integer n, sig real*8 body(n,2), x, y real*8 eps parameter(eps=1.d-3) if ((x.ge.body(1,1)-sig*eps).and.(x.le.body(2,1)+sig*eps).and. : (y.ge.body(1,2)-sig*eps).and.(y.le.body(3,2)+sig*eps)) then inside_polygon = .true. else inside_polygon = .false. endif return end c----------------------------------------------------------------------- c c Prusecik roviny (urcene bodem a 2 vektory: rp = r0 + Ap*a + Bp*b) c a primky (urcene bodem a vektorem: rp = r + Cp*c). Jeslize jsou c rovnobezne, lflag na .false. c subroutine intersect_plane_line(r0,a,b,r,c,Ap,Bp,Cp,lflag) implicit none real*8 r0(3),a(3),b(3),r(3),c(3),Ap,Bp,Cp logical lflag integer i real*8 a_gem(3,3), b_gem(3) logical sing do i = 1, 3 a_gem(1,i) = a(i) a_gem(2,i) = b(i) a_gem(3,i) = -c(i) b_gem(i) = r(i) - r0(i) enddo call gem(a_gem, b_gem, 3, 3, sing) if (.not.sing) then Ap = b_gem(1) Bp = b_gem(2) Cp = b_gem(3) lflag = .true. else lflag = .false. endif return end c----------------------------------------------------------------------- c c Otoceni 3D vektoru do roviny ciselniku (pomoci 2 matic rotace). c subroutine rotat0(x,az,ht) implicit none real*8 x(3),az,ht real*8 x0(3) call rotat2(x,x0,ht) call rotat3(x0,x,az) return end c----------------------------------------------------------------------- c c Vypocet A, h stredniho slunce (pro dane ekliptikalni souradnice l, b) c a pozorovaci stanoviste. c subroutine stredni_slunce_ah(l,b,t,lambda,phi,a,h) implicit none real*8 l,b,t,lambda,phi,a,h real*8 t0,x(3),x0(3),s0,ss,s,epsilon,alpha real*8 pi,degrad parameter(pi = 3.1415926535d0, degrad = pi/180.0d0) c functions real*8 nula2pi,eps_earth x0(1) = -cos(l)*cos(b) x0(2) = -sin(l)*cos(b) x0(3) = -sin(b) t0 = 0.d0 epsilon = eps_earth(t0) call rotat1(x0,x,epsilon) alpha = atan2(x(2),x(1)) c delta = asin(x(3)) c hvezdny cas se nepocita z JD, ale odvozuje z delky!!! <- toto bylo chybne!!! cc call lst(jd,lambda/degrad,s0,ss,s) c s0 = l/pi*12.d0 c spravne se ma odvozovat z rektascenze stredniho slunce, aby pro t=12 bylo na jihu!!! s0=alpha/pi*12.d0+12.d0 c zde 1 h slunecniho casu odpovida 15 stupnum v hodinovem uhlu!!! cc ss=s0+1.002737909350795d0*t ss=s0+t s=ss+lambda/degrad/15.d0 call rotat3(x,x0,s*pi/12.d0) call rotat2(x0,x,pi/2.d0-phi) a = nula2pi(-atan2(x(2),x(1))) h = asin(x(3)) return end c----------------------------------------------------------------------- c c Vypocet casu vychodu a zapadu pro stredniho slunce (o danem l, b) c a pozorovaci stanoviste (phi). c subroutine t_vychod_zapad(l,b,lambda,phi,t_vychod,t_zapad) implicit none real*8 l,b,lambda,phi,t_vychod,t_zapad integer i real*8 t,t0,x(3),x0(3),epsilon,z,s,alpha,sindelta,cosphicosdelta, : arccos,s0,s_vz(2),t_vz(2) real*8 pi,degrad,eps parameter(pi = 3.1415926535d0, degrad = pi/180.d0, eps = 1.d-38) c functions real*8 nula2pi,eps_earth,ut1 x0(1) = cos(l)*cos(b) x0(2) = sin(l)*cos(b) x0(3) = sin(b) t0 = 0 epsilon = eps_earth(t0) call rotat1(x0,x,-epsilon) alpha = atan2(x(2), x(1)) sindelta = x(3) cosphicosdelta = cos(phi) * sqrt(1.d0 - sindelta**2) s0 = l/pi*12.d0-12.d0 t_vychod = -1.d0 t_zapad = -1.d0 z = 90.d0*degrad if (abs(cosphicosdelta).gt.eps) then arccos = (cos(z) - sin(phi)*sindelta)/cosphicosdelta if (abs(arccos).le.1.d0) then t = acos(arccos) do i = 1, 2 s_vz(i) = alpha + (2*i-3)*t s_vz(i) = s_vz(i)/pi*12.d0 c zde se opet NEpocita korektne cas UT1, ale vezme se primo rozdil hvezdnych casu c the time UT1 is NOT calculated correctly here, this is only a difference of sidereal times c t_vz(i) = ut1(s_vz(i),t0,???,epsilon,lambda) t_vz(i) = s_vz(i) - s0 t_vz(i) = nula2pi(t_vz(i)/12.d0*pi)/pi*12.d0 enddo t_vychod = t_vz(1) t_zapad = t_vz(2) endif endif return end c----------------------------------------------------------------------- c c Kontrola, zda alespon NEKDY behem roku slunce sviti nad ciselnikem c pro rysky pocitane od vychodu nebo od zapadu c Check, if at least SOMETIMES during the year, the Sun is above the dial c for hour counted from sunrise or from sunset logical function t_vychod_zapad_chk(t0,lambda,phi,r0,a0,b0, : r_nodus,t_horiz,v_or_z) implicit none real*8 pi, pi2, deg, rad parameter(pi = 3.1415926535d0, pi2 = pi/2.d0, deg = pi/180.d0, : rad = 180.d0/pi) real*8 t0,lambda,phi,r0(3),a0(3),b0(3),r_nodus(3) integer t_horiz,v_or_z real*8 b,l,l0,t_vychod,t_zapad,t,a,h,c0(3),Ap,Bp,Cp,eps integer n,j logical lhour,lflag lhour = .false. l0 = 90.d0*deg eps = 1.d-3 b = 0. n = 180 j = 0 10 continue l = l0 + j*180.0d0*deg/n call t_vychod_zapad(l,b,lambda,phi,t_vychod,t_zapad) if (v_or_z.eq.1) then t = t0 + t_vychod else t = t0 + t_zapad endif call stredni_slunce_ah(l,b,t,0.d0,phi,a,h) c0(1) = cos(a)*cos(h) c0(2) = sin(a)*cos(h) c0(3) = sin(h) call intersect_plane_line(r0,a0,b0,r_nodus,c0, : Ap,Bp,Cp,lflag) c velke kladne Cp je zde kvuli rysce pro 6 h na J hodinach c large positive Cp is here due to the line for 6th hour on a southern dial if (((Cp.lt.0d0).or.(Cp.gt.1.d8)).and. : ((t_horiz.eq.1).or.(h.gt.-eps))) lhour = .true. j = j+1 if ((j.le.n).and.(.not.lhour)) goto 10 t_vychod_zapad_chk = lhour return end c----------------------------------------------------------------------- c c Interpolace mezi 2 realnymi cisly. c real*8 function linterp1(x,y,t,iflag) implicit none integer iflag real*8 x,y,t if (iflag.eq.0) then linterp1 = t*x + (1.d0-t)*y else linterp1 = (1.d0-t)*x + t*y endif return end c---------------------------------------------------------------------- c c Znamenko cisla (0 ma 1). c real*8 function plusminus(x) implicit none real*8 x if (x.ge.0.d0) then plusminus=1.d0 else plusminus=-1.d0 endif return end c---------------------------------------------------------------------- c c Ekliptikalni delka z delky dne (b = 0). c subroutine delka_dne_l(d,phi,l,iflag) implicit none real*8 d,phi,l integer iflag real*8 pi parameter(pi = 3.1415926535d0) real*8 t0, epsilon, de, sinl c functions real*8 eps_earth t0 = 0.d0 epsilon = eps_earth(t0) de = atan(-cos((d/12.d0*pi)/2.d0) / tan(phi)) sinl = sin(de)/sin(epsilon) if (abs(sinl).le.1.d0) then iflag = 0 l = asin(sinl) else iflag = 1 endif return end c---------------------------------------------------------------------- c c Delka dne z ekliptikalni delky (a sirky). c real*8 function l_delka_dne(phi,l,b) implicit none real*8 phi,l,b real*8 pi,deg parameter(pi = 3.1415926535d0, deg = pi/180.d0) real*8 t0,epsilon,sinde,cosde,z,eps,cosphicosde,cosd2,de c functions real*8 eps_earth t0 = 0.d0 epsilon = eps_earth(t0) sinde = sin(b)*cos(epsilon) + cos(b)*sin(l)*sin(epsilon) cosde = sqrt(1.d0-sinde**2) cosphicosde = cos(phi)*cosde eps = 1.d-16 de = asin(sinde) if (((phi.gt.0.d0).and.(de.gt.(90.d0*deg-phi))).or. : ((phi.lt.0.d0).and.(de.lt.(-90.d0*deg-phi)))) then l_delka_dne = 24.d0 else l_delka_dne = 0.d0 endif if (abs(cosphicosde).lt.eps) then return endif z = 90.d0*deg cosd2 = (cos(z)-sin(phi)*sinde)/cosphicosde if (abs(cosd2).gt.1.d0) then return endif l_delka_dne = 2.d0*acos(cosd2) / pi*12.d0 return end c---------------------------------------------------------------------- c c Vynechani znaku v retezci. c subroutine str_trim(s,c) implicit none character*(*) s character c integer i,j c functions integer length j = 0 do i = 1, length(s) if (s(i:i).ne.c) then j = j + 1 s(j:j) = s(i:i) endif enddo s = s(1:j) return end c---------------------------------------------------------------------- c c Zamena znaku v retezci za jiny. c subroutine str_replace(s,c1,c2) implicit none character*(*) s character c1,c2 integer i c functions integer length do i = 1, length(s) if (s(i:i).eq.c1) then s(i:i) = c2 endif enddo return end c********************************************************************** subroutine hhms24(h2,h,m,s) c********************************************************************** c c conversion of decimal hours to hours (in 0-24 interval), minutes and seconds c double precision h2,h,m,s,y double precision frac y=frac(h2/24.d0)*(24.d0+1.d-3) if (y.lt.0d0) y=y+pi2 h=int(y) m=int((y-h)*60) s=(y-h)*3600-m*60 return end c=======================================================================