	Program CYKLE
c 23 VIII 2008 poprawiono: date(today), real*8 ce,se w EARTH; 
c 2 IX nowe procedury Ksiezy i Earth-> Sun
c 3 X 2008 uaktualniono dtuts: do 2012 wg AA 2009
c v.2.1, 16 VI 2009 dodana funkcja narg i sprawdzanie liczby wprowadzanych danych 
c v. 2.2, 5 IX 2012 - dodane objasnienia na ekranie po starcie,
c                     ograniczenie liczby lat do 100.

	character*10 input,plik*30,file*30,znak(13)*1,z(0:1)*1,today*8
     &, rec*78,mies(12)*4,tydzien(7)*12,znakiZ(12)*10,wersja*3
c	IMPLICIT REAL*8 (D)
	real*8 dj0,djm,dj,dji(620),hpdeg,synodm	!
	real dl(620),hm(14),hp(14)
	integer r(14),m(14),d(14),ty(14),g(14),mi(14),zz(14)
	data z/'M','F'/
	data mies/' I',' II ','III ',' IV ',' V ',' VI ',
     &'VII ','VIII',' IX ',' X ',' XI ','XII'/,
     & znakiZ/'Baran','Byk','Bliznieta','Rak','Lew','Panna',
     & 'Waga','Skorpion','Strzelec','Koziorozec','Wodnik','Ryby'/,
     & tydzien/'Poniedzialek','Wtorek','Sroda','Czwartek','Piatek'
     &,'Sobota','Niedziela'/
     &,Wersja/'2.2'/

C Get current system date and time
c	CALL getdat(iy,im,id)
c*	call time(czas)	!  hh:mm:ss:ss
	call date(today)	!  MM/DD/YY
c The following is to amend a Y2k problem in F77L v.2.2
	if(ichar(today(7:7)).gt.57) today(7:7)=char(ichar(today(7:7))-10)
	read(today,'(i2,2(1x,i2))') imt,idt,iyt
	iyt=iyt+2000

	print'(32x,a//a/a//a,a/)',	'Program CYKLE',	!/31x,a,i5,2(1h.,i2.2)
c     &  '(v. 2009.06.16)',
     * '    Sluzy do obliczania dat fazy Ksiezyca - takiej jaka byla w p
     *odanym dniu (np. urodzenia sie kobiety). Wyniki zawieraja tez okre
     *slenie rodzaju znaku zodiaku:  M = meski (Baran, Bliznieta, ...)  
     *lub  F = zenski (Byk, Rak, ...).',
     *'    Obliczone dane stanowia kalendarz dni plodnych i plci dziecka
     * w tzw. CYKLU   KSIEZYCOWYM opisanym w "Nieznanym Swiecie"  (Nr 7-
     *8/1994, str. 71-73).          Pelny opis na stronie: '//
     &'www.astro.uni.torun.pl/~kb/Efemerydy/Cykle.htm'   
     &,'    <Wersja '//Wersja//' kompilowana dnia 2012.09.05',	!iyt,imt,idt,
     &'; K.M. Borkowski, CA UMK, Torun>'

c Czas [h] przejscia Ksiezyca kata x [deg] wynosi x*hpdeg (54,643164 h w znaku)
	hpdeg=27.321582d0*24d0/360	! Miesiac zwrotnikowy w godzinach/360 
	synodm = 29.530589d0 	! miesiac synodyczny
	plik='Cykle.txt'

	ir=1982
	im=04
	id=20
	igm=1341
	istr=1
	istrUrodz=2
	LetniCzas=1

30	print*,'    Wpisz wymagane dane; najpierw data i godz. urodzenia '
     &//'kobiety.  Nalezy podac kolejno: rok, miesiac, dzien, godzine i
     & Nr strefy czasowej urodzenia w formacie'
c	print*,
     &//'yyyy,MM,dd,hhmm,Nr.  Nr strefy'
     &//' to: 0 - UTC (GMT), 1 - nasz czas zimowy (CSE),   2 - letni '
     &//'(CWE), ... (dla USA i Kanady Nr jest ujemny, -4 do -8).      P'
     &//'rzyklad pierwszego wpisu: "1982 4 20 1341 2" oznacza 1982 r., k
     &wiecien 21, godz. 13:41  strefy +2  (tj. CWE). Kolejne wpisy dotyc
     &za wynikow: zakresu lat (dwie liczby), strefy i czasu letniego (dw
     &ie cyfry; tu zwykle mozemy nic nie pisac, samo Enter)i nazwa pliku
     & z wynikami (ciag znakow bez spacji i przecinkow).'
	print*
     &,'  UWAGA: po wpisaniu nowych danych program powtorzy zapytanie (w
     &yswietlajac nowe  dane wewnatrz nawiasow kwadratowych) i wtedy, je
     &sli dane sa OK, akceptujemy je  nacisnieciem klawisza Enter.  
     &  Teraz beda cztery zapytania o dane.'
	print*
31	print'(a,i5,2i3.2,i5.3,i3,a)','  Rok, mies., dzien, hhmm, Nr'//
     & ' [ENTER =',ir,im,id,igm,istrUrodz,']: '
	read(*,'(a)') rec
		if(rec.ne.' ') then
	n=narg(rec,last)
		if(n.lt.5) then
	  rec(last+1:min0(78,last+6))=',,,,,'
	  print*,' Wpisano za malo danych: "'//rec(1:last)
     &           //'" - sprawdz czy jest OK!'
		endif
	read(rec,*) ir,im,id,igm,istrUrodz
	go to 31
		endif
	min=mod(igm,100)
	ig=igm/100
		if(ir.lt.0) then ! nigdy
	print'(a,i4,a)','  Rok [ENTER=',ir,']: '
	read(*,'(a)') input
	if(input.ne.' ') read(input,'(i4)') ir
	print'(a,i2,a)','& Miesiac [ENTER=',im,']: '
	read(*,'(a)') input
	if(input.ne.' ') read(input,'(i2)') im
	print'(a,i2,a)','& Dzien [ENTER=',id,']: '
	read(*,'(a)') input
	if(input.ne.' ') read(input,'(i2)') id
	print'(a,i2.2,i2.2,a)','& Godz.,min. [ENTER=',ig,min,']: '
	read(*,'(a)') input
c	if(input(1:2).ne.' ') read(input(1:2),'(i2)') ig
c	if(input(4:5).ne.' ') read(input(4:5),'(i2)') min
	if(input.ne.' ') read(input,'(2i2)') ig,min
	print'(a/a,i2,a)','& Podaj Nr strefy czasowej <0=GMT, 1=CSE '
     *//'(zimowy w Polsce), 2=CWE (letni),  ','    3=moskiewski itd.> [E
     &NTER=',istrUrodz,']: '
	read(*,'(a)') input
	if(input.ne.' ') read(input,'(i2)') istrUrodz !,LetniCzas
	print'(a,i5,4i3.2,i3,a)','  Wczytano dane:',ir,im,id,ig,min,
     *istr,'.  Czy nie ma bledu [ENTER = Nie ma]? '
	read(*,'(a)') input
	if(input.ne.' ') go to 30
		endif

	etut=dtuts(ir+0.5)/(3600*24.)
	dj0=JD(ir,im,id,0)-0.5d0
	dj0=dj0+(ig-istrUrodz+min/60d0)/24d0+etut
c Faza Ksiezyca i znak zodiaku na moment poczatkowy
	CALL Sun(DJ0,RA,DEC,DS,Radius,EQT,gst)
	CALL KSIEZY(DJ0,RA,DEC,DK,SK,Radius)
	FAZA0=AMOD((DK-DS)/6.283187+1.,1.)
c	dlugK=dk/6.283187*360

	iy1=2008	!max0(ir+15,1994)
	iy2=iy1+2
40	print'(a,2i5,a)',
     *'  Na jakie lata (od do) liczyc kalendarze [ENTER =',iy1,iy2,']? '
	read(*,'(a)') rec
		if(rec.ne.' ') then
	n=narg(rec,last)
c	print*,n,last
			if(n.lt.2) then
	  rec(last+1:min0(78,last+6))=',,,'
	  print*,' Wpisano za malo danych: "'//rec(1:last)
     &           //'" - sprawdz czy jest OK!'
			endif
	read(rec,*) iy1,iy2
	if(iy2-iy1.gt.99) then
	  print*,'   Uwaga: program ograniczyl liczbe lat do 100!'
	  iy2=iy1+99
	endif
	go to 40
		endif

42	print'(a/a,i3,i2,a)','  Nr strefy czasowej wynikow oraz '	!&
     *//'indeks czasu letniego <0 - tylko zimowy, ',
     &'  1 - uwzglednic miedzy ostatnimi niedzielami III i X> [ENTER =',
     &istr,LetniCzas,']: '
	read(*,'(a)') rec
		if(rec.ne.' ') then
	n=narg(rec,last)
			if(n.lt.2) then
	  rec(last+1:min0(78,last+6))=',,,'
	  print*,' Wpisano za malo danych: "'//rec(1:last)
     &           //'" - sprawdz czy jest OK!'
			endif
	read(rec,*) istr,LetniCzas
	go to 42
		endif



c	print*,' OK, Go !!!!'
10	do 11 lplik=30,1,-1
	if(plik(lplik:lplik).ne.' ') go to 12
11	continue
c	lplik=len(plik)
12	print'(a,a,a)','  Nazwa pliku na wyniki [ENTER = ',
     *plik(1:lplik),']: '
	read(*,'(a)') file
	if(file.eq.' ') go to 20
	plik=file
	go to 10
20	open(1, file=plik,status='unknown')

	print'(a,f8.5,a,f8.2,a/)','  Faza poczatkowa [liczona od 0 do 1]'
     & //' wynosi',faza0,'  czyli',faza0*360,' stopni.'

	write(1,'(/14x,a,i5,1h-,i4/a,i5,2i3.2,a,i3.2,1h:,i2.2,a,
     *i3/a/a,f6.1,a,i3,a,i3,a/a/a/)')
     * ' KALENDARZ CYKLU KSIEZYCOWEGO na lata',iy1,iy2,
     *' dla kobiety urodzonej w dniu',ir,im,id,' o godz.',ig,min,
     *' czasu strefy Nr',istrUrodz,
     *' [0=GMT, 1=CSE (zimowy w Polsce), 2=CWE (letni w Polsce) itd.]. F
     &aza Ksiezyca' ,
     *' w tym momencie wynosi',faza0*360,' stopni (ok.',
     *nint(faza0*synodm),' dni po albo',nint((1-faza0)*synodm),' dni prz 
     &ed nowiem).'
	write(1,'(a/a/a/a))')
     * ' W kazdym slupku podano datê: miesiac, nizej - dzien, nizej - go  
     *dzine i',' minuty (hh:mm), w czwartym wierszu - plec'//
     *' (M - chlopiec, F - dziewczynka)',' oraz znak zodiaku,'//
     &' a w piatym - liczbe godzin od momentu wejscia Ksiezyca ',
     &' w ten znak i do wyjscia z niego (w sumie, przed zaokragleniem,'
     &//' 54.6 godzin).'
		if(LetniCzas.eq.0) then
	write(1,'(23x,a,i3,a)') 'Daty i godziny podano w czasie strefy'
     & ,istr,'.'
		else
	write(1,'(a,i3,a,i3,a/)') ' Daty i godziny podano w czasie strefy'
     &,istr,' (czas zimowy) i',istr+1,' (letni).'
		endif
	dj=JD(iy1,1,26,0)
	djm=JD(iy2,12,31,0)+1.

	do 55 nf=1,620
c	if(mod(nf,13).eq.1) write(*,'(3x,''Rok'',i5)') iy1+nf/13
50	call czasF(DJ,FAZA0,dk)
	if(dj.gt.djm) go to 58
	dji(nf)=dj
	dl(nf)=dk*180/3.141593
55	dj=dj+31

58	i=0
c	print*,' Ilosc dat',nf
	do 100 iy=iy1,iy2
		if(LetniCzas.ne.0) then
	jdMar=JD(iy,3,31,0)
	ntydzien=mod(jdMar,7)	! 0 - poniedz., 6 - niedziela
	if(ntydzien.lt.6) jdMar=jdMar-(ntydzien+1)	! niedziela konca marca
	call DATA(JDMar,Lrok,Miesiac,NMar,0)
	jdOct=JD(iy,10,31,0)
	djMar=jdMar-.5d0 + 1d0/24	! 1:00 UTC (= 2:00 CET)
	ntydzien=mod(jdOct,7)
	if(ntydzien.lt.6) jdOct=jdOct-(ntydzien+1)	! ostatnia niedz. pazdz.
	call DATA(JDOct,Lrok,Miesiac,NOct,0)
	djOct=jdOct-.5d0 + 1d0/24	! 1:00 UTC (= 2:00 CET)
c	 print*,iy,' Mar:',NMar,' Oct',NOct
		endif
	etut=dtuts(iy+0.5)/(3600*24.)
	nm=0
60	nm=nm+1
62	dl(nm)=dl(i+nm)
	dj=dji(i+nm)+istr/24d0-etut
	  if(LetniCzas.ne.0.and.dj.gt.djMar.and.dj.lt.djOct ) dj=dj+1d0/24
	j=dj+0.5d0
	call DATA(j,r(nm),M(nm),d(nm),0)
	if(r(nm).gt.iy.or.i+nm.eq.nf) go to 65
	if(r(nm).lt.iy) i=i+1
	if(r(nm).lt.iy) go to 62
	godz=dmod(dj-.5d0,1d0)*24
	mi(nm)=amod(godz,1.)*60+.5
	ty(nm)=mod(j,7)+1
	g(nm)=godz
	angz=mod(dl(nm)+360.,30.) 	! kat pozycyjny Ksiezyca w znaku
	hm(nm) = angz*hpdeg
	hp(nm) = (30-angz)*hpdeg
	zZ(nm)=mod(int(dl(nm)/30.+12),12)+1
	iz=mod(zZ(nm)-1,2)
	znak(nm)=z(iz)
	go to 60
65	nm=nm-1
	i=i+nm
	write(*,'(a,i5)') '       Plec i daty dla roku',iy
c	write(*,'(1h ,1hl,i4,12(i6))') (nint(dl(k)),k=1,nm)
	write(*,'(13(1x,a1,2i2))') (znak(k),m(k),d(k),k=1,nm)
c	write(*,'(13(i3.2,1h:,i2.2))') (g(k),mi(k),k=1,nm)
		if(LetniCzas.eq.0) then
	write(1,'(/25x,a,i5)') 'Kalendarz dni plodnych w roku',iy
		else
	write(1,'(/7x,a,i5,a,i3,a,i3,a)') 'Kalendarz dni plodnych w roku',
     & iy,' (czas letni od',NMar,' III do',NOct,' X)'
		endif
c	write(1,'(1h ,1hl,i4,12(i6))') (nint(dl(k)),k=1,nm)
	write(1,'(13(2x,a4))') (mies(m(k)),k=1,nm)
	write(1,'(13(i3,1x,a2))') (d(k),tydzien(ty(k)),k=1,nm)
	write(1,'(13(i3.1,1h:,i2.2))') (g(k),mi(k),k=1,nm)
	write(1,'(13(1x,a1,1x,a3))') (znak(k),znakiZ(zZ(k)),k=1,nm)

	write(rec,'(13(i3,i3))') 
     & (nint(hm(k)),nint(hp(k)),k=1,nm)
	do 66 k=1,nm
	l=(k-1)*6+2
	if(hm(k).lt.0.95) write(rec(l:l+1),'(f2.1)') hm(k)
66	if(hp(k).lt.0.95) write(rec(l+3:l+4),'(f2.1)') hp(k)
	write(1,'(a)') rec
100	continue
	write(1,'(1h ,77(1h-))')
c	print*,'  Konczymy obliczenia [ENTER = tak] ? '
c	read(*,'(a)') input
c	if(input.ne.' ') write(1,*)
c	if(input.ne.' ') go to 30

	write(1,111) Wersja,iyt,imt,idt
111	format(/17x,'Obliczenia wykonano programem CYKLE ',a,
     *' (KMB) dnia',i5,2(1h.,i2.2))
	close(1)
	print*
	print*,'         Wyniki znajduja sie w pliku "',plik(1:lplik),'"'
	end


	function narg(input,last)
c Liczy ile jest argumentow oddzielonych przecinkiem lub spacja oraz
c zwraca miejsce ostaniego znaku nie bedacego spacja ani przecinkiem
	character*78 input
	n0=78
	narg=0
	i=n0+1
10	i=i-1
	if(i.eq.0) go to 1000
	if(input(i:i).ne.' '.and.input(i:i).ne.',') then
15	  narg=narg+1
	if(narg.eq.1) last=i
20	  i=i-1
	  if(i.eq.0) go to 1000
	  if(input(i:i).ne.' '.and.input(i:i).ne.',') go to 20
	endif
	go to 10
1000	end

	SUBROUTINE czasF(DJD,faza0,dk)
	DOUBLE PRECISION DJD
	DD=0.
1	CALL Sun(DJD-DD,RA,DEC,DS,R,EQT,gst)
	CALL KSIEZY(DJD-DD,RA,DEC,DK,SK,R)
	FAZ=AMOD((DK-DS)/6.283187+1.,1.)
	FAZ=faza0-FAZ
	IF(DD.EQ.0..AND.FAZ.GT.0.) FAZ=FAZ-1.
	IF(DD.NE.0..AND.ABS(FAZ).GT..5) FAZ=FAZ-SIGN(1.,FAZ)
	DD=DD-FAZ*29.530589d0
	IF(ABS(FAZ).GT..000025) GO TO 1
C  O.000025 FAZ'a odpowiada ok. 1 min. DD
	djd=djd-dd
	END


	subroutine ksiezy(DJ,RA,DEC,DLE,B,RK)

c   Procedura oblicza wspolrzedne [w rad] ekliptyczne (DLE i B) i
c  rownikowe (RA i DEC) oraz odleglosc ksiezyca od ziemi  [w promie-
c  niach ziemi] na moment okreslony przez dzien julianski (DJ).
c  Jest to uproszczona wersja MOON (blad 10"[dzis] - 20"[1500BC])
c  TIDACC = -26"/cy/cy mozna zmieniac (patrz -->PARAMETER)
	implicit real*8 (a-h,l,o-z)
	dimension ap1(4),al(5),bD(4),pl(4),bl(4),bf(4),eps(4)
	real*4 dle,b,ra,dec,rk,av(72),au(2,24),au3(25:32),ar(52)
	integer*2 iv(4,72),iu(2,24),iu3(25:32)
	parameter (rad=1/206264.8062470964d0,tidacc=-26.d0)
      Data  iv,iu,iu3 /0,0,1,0,2,0,-1,0,0,0,2,0,0,1,0,0,0,0,0,2,
     * 2, 0,-2, 0, 2,-1,-1, 0, 2, 0, 1, 0, 2,-1, 0, 0, 0, 1,-1, 0,
     * 1, 0, 0, 0, 0, 1, 1, 0, 2, 0, 0,-2, 0, 0, 1,-2, 4, 0,-1, 0,
     * 0, 0, 3, 0, 4, 0,-2, 0, 2, 1,-1, 0, 2, 1, 0, 0, 1, 1, 0, 0,
     * 2,-1, 1, 0, 2, 0, 2, 0, 4, 0, 0, 0, 1, 0, 1, 0, 2, 0, 0, 0,
     * 2, 0,-3, 0, 0, 1,-2, 0, 2,-1,-2, 0, 2,-2, 0, 0, 0, 1, 2, 0,
     * 4,-1,-1, 0, 2, 1, 1, 0, 1, 0,-1, 0, 2,-2,-1, 0, 0, 0, 2,-2,
     * 2, 0, 1,-2, 3, 0,-1, 0, 2, 2,-1, 0, 0, 2,-1, 0, 4,-1,-2, 0,
     * 1, 0,-2, 0, 4,-1, 0, 0, 4, 0, 1, 0, 0, 2, 1, 0, 0, 0, 4, 0,
     * 0, 2, 0, 0, 2,-1, 2, 0, 2,-2, 1, 0, 2,-1, 0,-2, 1,-1, 0, 0,
     * 3, 0, 0, 0, 0, 0, 1, 2, 2, 0,-1, 2, 2, 0, 0, 2, 0, 0, 2, 2,
     * 2, 1,-2, 0, 2, 1, 0,-2, 1, 1, 1, 0, 3, 0,-2, 0, 4, 0,-3, 0,
     * 1, 1,-1, 0, 2, 0, 3, 0, 2, 0, 1, 2, 2, 0,-4, 0, 0, 1,-3, 0,
     * 4, 1,-1, 0, 1, 0, 2, 0, 1, 0, 0,-2, 6, 0,-2, 0, 2, 0,-2,-2,
     * 0, 1, 3, 0, 2, 0,-2, 2,
     * 1,-1, 1,-1, 1,-1, 1,-1,-1,1,-1,1, 1,-1,-1,1,-1,1,-1,1, 1,-1,
     * 1,-1, 3,-1,-1,5,-1,1, 1,-1, 1,-1, 1,-1,-1,1,-1,1,-1,1,-1,1,
     *-1,1, 1,-1,        -1,-1,-1,-1,-1, 1,-1,-1/

      Data  av,au,au3,ar/22639.586,4586.44,769.03,0.,-411.60,211.66,
     *    0. ,191.96, 0.,0.,-124.99,0.,  55.18,39.53,38.43,36.12,
     *   30.77,-28.40,-24.36,17.95,14.53,14.38,13.9,-8.45,2369.91,
     *   13.19,-9.68,8.61,8.05,-7.63, 4.37,-2.91,-18.58,7.37,-1.37,
     *   -6.38,-3.21,-2.52,-2.57, 2.73,-1.75, 1.87,1.98,-1.16, 1.93,
     *    -7.45, 1.18, 0.75, 2.15,-0.56,.4,-45.10,-9.37,-5.74,-4.00,
     *    2.49,-1.44, 1.26,-1.22, 1.19, 1.08, 1.06,-0.99,0.95,-0.67,
     *    -0.64,-0.58,-0.58, 0.57,-0.56,-0.55,-0.54,
c au
     * 1010.17,999.69,199.48,166.57,61.91,31.76,-6.46,-4.84,
     *18461.24,-6.30,15.57,-1.62, 8.87, 7.43,33.36,15.12,
     *   29.58, 7.96,-6.73,-5.63,-5.37,-4.81,-5.31,-5.08,
     *  117.26, 2.19, 2.80,-1.02, 6.58, 3.00, 3.98, 1.58,
     *    2.41, 0.63,-1.32,-0.79,-12.09,-1.26,0.80, 0.80,
     *    1.77, 1.13, 2.15, 1.52, 3.67, 1.19,-0.67,-0.59,
c au3
     *  623.65, 1.52,-0.79, 0.65, 1.09,-0.64, 0.60,-0.82,
c ar
     *-20905.36, -3699.11,  -569.93,   0.  ,   -3.15,246.16,
     *     .0  ,  -170.73,     0.  ,   0.  ,  108.74,  0.  , 10.32,
     *    79.66,   -34.78,   -23.21, -21.64,   24.21, 30.82,-16.67,
     *   -12.83,   -10.44,   -11.65,   6.32,-2955.97, 14.40, -7.00,
     *    10.06,-9.88, 5.75,-3.96, 2.62,-8.38,-4.95,
     *    -4.42, 4.13, 3.26, 2.35,-2.12,-1.90,-1.74,
     *    -1.57,-1.42, 1.17,-1.12, 1.07,-0.85,-0.66,
     *     0.66, 0.50,-1.42, 0.00/
c zgodne z LE51 (Tab 8,AA 190,342); do AL (Lk) dodano 0.0351"*t**2
        data al/0.d0,1732564372.29232d0,-5.9462d0,0.6642d-2,-0.5522d-4/
        data bD,pl,bl,bf/1602961600.8820d0,-7.0731d0,0.6556d-2,
     1  -0.3184d-4,  129596581.0996d0,-0.5529d0,0.147d-3,0.d0,
     2  1717915922.8022d0,31.1665d0,0.51612d-1,-0.24470d-3,
     3  1739527262.7141d0,-13.4508d0,-0.1060d-2,0.417d-5/
        data ap1/210664136.43355d0,
     1  129597742.2758d0,68905077.59284d0,10925660.42861d0/
     *,eps/84381.448d0,-46.815d0,-.00059d0,.001813d0/
	dms(id,im,s)=(id*60+im)*60+s

	t=(dj-2451545.D0)/36525.D0
	dk=(tidacc+26.3046d0)*(.5+1.60816d-5*t)
	Lk=(dms(218,18,59.95571D0+.24177d0)+
     *  (((al(5)*t+al(4))*t+al(3)+dk)*t+al(2))*t)*rad
        D=(dms(297,51,0.73512D0 -.07293d0)+
     * (((bD(4)*t+bD(3))*t+bD(2)+dk)*t+bD(1))*t)*rad
        lp=(dms(357,31,44.79306D0+.03845)+
     * (((pl(4)*t+pl(3))*t+pl(2))*t+pl(1))*t)*rad
        l=(dms(134,57,48.28096D0 -.097)+
     * (((bl(4)*t+bl(3))*t+bl(2)+dk)*t+bl(1))*t)*rad
        F=(dms(93,16,19.55755D0 -.05596d0)+
     * (((bF(4)*t+bF(3))*t+bF(2)+dk)*t+bF(1))*t)*rad

        Ve=(dms(181,58,47.28305D0)+ap1(1)*t)*rad
        Te=(dms(100,27,59.22059D0)+ap1(2)*t)*rad
	aMa=(dms(355,25,59.78866d0)+ap1(3)*t)*rad
        aJu=(dms(34,21,5.34212D0)+ap1(4)*t)*rad
	av(4)=-666.4171+t*(1.6768+t*.004787)
	av(7)=205.43582-t*(.51642+t*.0015)
	av(9)=164.72851-t*(.41383+t*.0012)
	av(10)=-147.32129+t*(.37115+t*.00108)
	av(12)=-109.38029+t*(.2756+t*.0008)
	ar(9)=-204.59+t*(.51395+t*.00149)
	ar(7)=-152.14+t*(.38245+t*.00111)
	ar(10)=-129.62+t*.32654
	ar(12)=104.76-t*.26396
	ar(4)=48.8883-t*.12302

	VeT=18*Ve-16*Te-l+.46326
	sV=sin(VeT)
	r=8.752*cos(D+D-l-F-F)+1.059*cos(D+D-l+Te+Te-aJu-aJu+.0021)
	u=-8.0451*sin(Lk)+1.5102*sin(D+Te+4.829)+1.2605*sV*cos(F)
     * -.87*sin(l)*cos(Lk)-.614*sin(Lk-F-F)*cos(l+l-F)
     * -.63*sin(l+l)*cos(D+D-F-F-F)
	v=sV*(14.24883+1.5274*cos(l))+(7.063+.9847*cos(Lk-F+l+l))*
     * sin(Lk-F)-1.14307*sin(D+D-l+Te+Te-aJu-aJu)  +.90114*sin(4*Te
     * -8*aMa+3*aJu+4.9914)    +(.603*cos(Ve-Te)-.82155)*sin(Ve-Te)
     * -.6437*sin(3*(Ve-Te)+D+D-l)   +.6388*sin(Te-aJu+.02145)
     * +.56341*sin(10*Ve-3*Te-l+5.81728)     +t*.25425*cos(VeT)

	do 3 i=33,52
	j=i+20
	df=iv(1,i)*D+iv(2,i)*lp+iv(3,i)*l+iv(4,i)*F
	v=v+av(i)*sin(df)+
     *    av(j)*sin(iv(1,j)*D+iv(2,j)*lp+iv(3,j)*l+iv(4,j)*F)
3	r=r+ar(i)*cos(df)

	do 2 i=25,32
	df=iv(1,i)*D+iv(2,i)*lp+iv(3,i)*l+iv(4,i)*F
	u=u+au3(i)*sin(iu3(i)*F+df)
	v=v+av(i)*sin(df)
2	r=r+ar(i)*cos(df)

	do 1 i=24,1,-1
	df=iv(1,i)*D+iv(2,i)*lp+iv(3,i)*l+iv(4,i)*F
	u=u+au(1,i)*sin(iu(1,i)*F+df)+au(2,i)*sin(iu(2,i)*F+df)
	v=v+av(i)*sin(df)
1	r=r+ar(i)*cos(df)

	rk=(r+385000.529d0)/6378.14d0
	b=u*rad
	dla=-(17.2+.0174*t)*sin(Lk-F)-1.3187*sin(Lk-D+Lk-D)-
     * .43*sin(F)*cos(Lk+Lk-F)
	dle=dmod(Lk+(v+dla)*rad,6.283185307178d0)
		 if(dle.lt..0) dle = dle+6.2831853
	ep=((eps(4)*t+eps(3))*t+eps(2))*t+eps(1)
	se=sin((ep+9.2025*cos(Lk-F)+.5736*cos(Lk-D+Lk-D))*rad)
	ce = sqrt(1. - se*se)
	sd = sin(dle)
	dec= asin(sd*se*cos(b) + sin(b)*ce)
	ra=dmod(datan2(sd*ce-tan(b)*se,dcos(dle+0d0))+6.2831853d0,
     *6.2831853d0)
	end


	subroutine sun(DJ,RA,DEC,DL,RS,EQT,GST)

C  PROCEDURA OBLICZA WSPOLRZEDNE ROWNIKOWE (RA I DEC), DLUGOSC 
C  EKLIPTYCZNA (DL) [rad] ORAZ ODLEGLOSC SLONCA OD ZIEMI (R [AU])
c  Algorytm: P.Bretagnon, J.L.Simon, J.Laskar,J.Hist.Astron.17(1986),39-50
c  blad dlug (dld) = 2"+.03"*(u*10)**2
	REAL*8 DJ,u,rsd,dld,f,a(38),n(38),l(38),r(18),pi2,t
	COMMON /SUNCOM/ XI0,ETA0,SE,CE,DLa /PLAN/ T

	data pi2/6.283185307179D0/
     *, a/6.24005d0,6.1969147d0,6.15378d0,4.315d0,5.198d0,2.846d0,
     * 1.423d0,8.63d0,3.193d0,1.223d0,2.75d0,9.944d0,4.5d0,2.84d0,
     * 1.92d0,4.27d0,1.89d0,5.98d0,4.533d0,0.061d0,2.828d0,4.654d0,
     * 3.229d0,4.374d0,4.345d0,3.44d0,4.24d0,5.96d0,.09d0,4.03d0,
     * 2.65d0,1.72d0,4.27d0,.93d0,2.21d0,3.59d0,4.97d0,5.69d0/

     *, n/3*.0d0,57533.85d0,777137.71d0,78604.2d0,39302.1d0,115067.7d0,
     * 15774.34d0,15773.85d0,52237.69d0,58849.26d0,55076.5d0,55075.7d0,
     * 54868.6d0,117906.3d0,109771.2d0,-55731.4d0,-33.93d0,-34.86d0,
     * 5296.67d0,5296.11d0,261.08d0,264.89d0,-3980.7d0,-7756.6d0,
     * -7752.8d0,-7961.4d0,25443.9d0,60697.8d0,207.8d0,2132.2d0,2132.8d0
     *,-8.d0,46941.1d0,-68.3d0,29424.6d0,157208.4d0/

     *, l/3*0.d0,350.d0,314.d0,-268.d0,234.d0,132.d0,129.d0,64.d0,78.d0,
     * -99.d0,72.d0,29.d0,24.d0,-32.d0,27.d0,21.d0,334.d0,158.d0,114.d0,
     * 93.d0,68.d0,37.d0,86.d0,38.d0,14.d0,28.d0,21.d0,20.d0,13.d0,
     * 27.d0,18.d0,12.d0,2*10.d0,13.d0,-10.d0/

     *, r/3*0.d0,-163.d0,309.d0,158.d0,-54.d0,-93.d0,-23.d0,-11.d0,
     * -33.d0,47.d0,-33.d0,-14.d0,-11.d0,24.d0,-19.d0,31.d0/

	T  = DJ-2451545d0
	u=(dj-2451545D0)/3652500.d0
	l(1)=((2885d0*u-25347d0)*u-84065d0)*u+334166d0
	l(2)=((194d0*u-309d0)*u-1755d0)*u+3489d0
	l(3)=51-38*u
	r(1)=-.499961*l(1)
	r(2)=-.4*l(2)
	r(3)=-.346*l(3)
	n(1)=(((-.0024*u-.0055)*u+.0007)*u-.02682)*u+62830.1955d0
	n(2)=2.*n(1)
	n(3)=3.*n(1)
	dld=(dmod(4.8950592d0+u*62833.1966661d0,pi2)+u*u*(((((.00063*u+
     * .00082)*u-.00088)*u-.011408)*u+.00035)*u+.052919))*1d7
	rsd=((80*u-120)*u-702d0)*u
	a7=dld*1d-7
	do 1 i=38,1,-1
	f=dmod(a(i)+n(i)*u,pi2)
	if(i.lt.19) rsd=rsd+r(i)*dcos(f)
1	dld=dld+l(i)*dsin(f)
c do nutacji: l(1) = omega, l(2) = 2Ls, dlk2 = 2Lk
	l(1)=dmod(2.1824d0+(.36244*u-3375.7)*u,pi2)
	l(2)=dmod(3.507d0+(.106*u+125666.39)*u,pi2)
	dlk2=dmod(2*(a7+a(5)+n(5)*u),pi2)
	dla=((-17.1996-1.742*u)*dsin(l(1))-1.3187*dsin(l(2))-
     -.2274*sin(dlk2)+.2062*dsin(2*l(1)))/206264.8
c  u gory nutacja w dlug.[rad];        nizej - nutacja w eps ["]
	de=9.2025*dcos(l(1))+.5736*dcos(l(2))+.0977*cos(dlk2)-
     -.0895*dcos(2*l(1))
	se=dsin((23+26/60d0+(((((((10.*u-39.)*u-250.)*u-51.4)*u+
     *    1999.3)*u-1.6)*u-4680.93)*u+21.448+de)/3600.)*pi2/360)
2	ce=sqrt(1-se*se)
c aberracja = -(993.651+17.*cos(f{i=1}))*1d-7 [rad]
	dld=dmod((dld-(993.651+17.*dcos(f)))*1d-7+dla+pi2,pi2)
	dl=dld
	rs=rsd*1d-7+1.0001399d0
	XI0 = Rs*COS(DLd)
	ETA0= Rs*SIN(DLd)
	RA  = ATAN2(CE*sin(dld),cos(dld))
	DEC = ASIN(SE*sin(dld))
	EQT = AMOD(a7 - RA,6.283185)
	 IF(ABS(EQT).GT.3.)  EQT = EQT - SIGN(6.2831853,EQT)
	 IF(RA.LT..0) RA = RA + 6.2831853
	GST = DMOD((DJ-2451545D0)*pi2+A7+DLA,pi2)
	 IF(GST.LT.0.) GST = GST + 6.2831853
	END

	function dtuts(y)
	dimension idt(32),dt(233)		! dt(233=2012-1780+1)
c wg: FR Stephenson, LV Morrison, 1984, Phil Trans R Soc Lond A313, 47
c iy 1630 - 1775 (1625 dodalem wg rownania; 1780 z nast.)
	Data idt/78,    85,72,62,54,48,43,37,32,26,21,16,12,10,9,9,9,
     * 10,10,11,11,11,12,12,13,13,14,15,16,16,17,    17/
c iy 1780 - 1859
	Data dt/2*16.9,17.,6*17.1,17.,16.9,16.7,16.5,16.2,15.9,15.6,15.2
     *,14.8,14.4,14.1,13.7,13.4,13.1,12.9,12.7,12.6,11*12.5,12.4,12.3,
     *12.2,12.,11.7,11.4,11.1,10.6,10.2,9.6,9.1,8.6,8.,7.5,7.,6.6,6.3,6.
     *,5.8,5.7,3*5.6,5.7,5.8,5.9,6.1,6.2,6.3,6.5,6.6,6.8,6.9,7.1,7.2,
     *7.3,7.4,7.5,7.6,2*7.7,2*7.8,
c iy 1860 - 1980
     *7.88,7.82,7.54,6.97,6.4,6.02,5.41,4.1,2.92,1.82,1.61,
     *.1,-1.02,-1.28,-2.69,-3.24,-3.64,-4.54,-4.71,-5.11,-5.4,-5.42,
     *-5.2,2*-5.46,-5.79,-5.63,-5.64,-5.8,-5.66,-5.87,-6.01,-6.19,
     *-6.64,-6.44,-6.47,-6.09,-5.76,-4.66,-3.74,-2.72,-1.54,-.02,1.24,
     *2.64,3.86,5.37,6.14,7.75,9.13,10.46,11.53,13.36,14.65,16.01,17.2,
     *18.24,19.06,20.25,20.95,21.16,22.25,22.41,23.03,23.49,23.62,
     *23.86,24.49,24.34,24.08,24.02,24.,23.87,23.95,23.86,23.93,23.73,
     *23.92,23.96,24.02,24.33,24.83,25.3,25.7,26.24,26.77,27.28,27.78,
     *28.25,28.71,29.15,29.57,29.97,30.36,30.72,31.07,31.35,31.68,
     *32.18,32.68,33.15,33.59,34.,34.47,35.03,35.73,36.54,37.43,38.29,
     *39.2,40.18,41.17,42.22,43.37,44.48,45.47,46.46,47.52,48.53,49.59,
     *50.54,
c iy 1981 - 2007 (uaktualniono 2008.10.03 wg. Astron.Almanac 2009)
c iy:  1981    82    83     84    85   86    87    88    89    90   91
     *51.38,52.17,52.96,53.79,54.34,54.87,55.32,55.82,56.30,56.86,57.57
cc     1992   93    94    95    96     97    98    99   2000   1
     *,58.31,59.12,59.98,60.78,61.63,62.29,62.97,63.47,63.83,64.09,
c                 (AA extrapolated to 2008 - 2012)
cc    2002   3     4     5     6    7     8    9  10     11   12
     *64.3,64.47,65.57,64.69,64.85,65.15,65.5,66.1,66.7 ,67., 68./
	iytab=2012	!wymiar tab. dt(ii): ii=iytab+1-1780

	iy=y
	f=y-iy
	if(iy.gt.iytab-1.or.iy.lt.1637) go to 2

	if(iy.lt.1780) go to 1
	dtuts=dt(iy-1779)*(1-f)+dt(iy-1778)*f
	return

1	i=(iy-1620)/5
cc	dtuts=idt(i)+(iy-i*5-1619.5)*(idt(i+1)-idt(i))/5.
	dtuts=idt(i)+(y-i*5-1620)*(idt(i+1)-idt(i))/5.
	return
2	t=(y-1800)/100.
c w zasadzie wstecz tylko do ok. 700 BC
	if(iy.lt.948) dtuts=(44.3*t+320)*t+1360.
	if(iy.ge.948) dtuts=25.5*t*t 	! 44 s dla ciaglosci; nizej

c musze odjac dla ciaglosci z danymi obserw. w XXI w.
		if(iy.ge.iytab) then
	dtutii=dt(iytab+1-1780)		! ostatnia wartosc w tabeli
	dtutform=25.5*((iytab-1800)/100.)**2 ! odpwiada jej ta wg wzoru
	dtuts= dtuts-(dtutform-dtutii)
		endif
	end

	FUNCTION JD(L,M,N,J1G0)

C   OBLICZA DZIEN JULIANSKI Z DATY KALENDARZA GREGORIANSKIEGO [J1G0=0]
C  LUB JULIANSKIEGO [J1G0=1]  (L = ROK, M = MIESIAC, N = DZIEN MIESIA-
C  CA). ALGORYTM DZIALA POPRAWNIE OD 1 MARCA -100100 R.  (OBU KALENDA-
C  RZY).

      JD=(L+(M-8)/6+100100)*1461/4+(153*MOD(M+9,12)+2)/5+N-34840408
      IF(J1G0.LE.0)          JD = JD-(L+100100+(M-8)/6)/100*3/4+752
	  END

	SUBROUTINE DATA(JD,L,M,N,J1G0)

C   OBLICZA DATE GREGORIANSKA LUB JULIANSKA  (L = ROK, M = MIESIAC
C  [STYCZEN: M=1], N = DZIEN MIESIACA)  Z DNIA JULIANSKIEGO (JD) W
C  POLUDNIE.  STOSOWALNOSC OD JD = -34839655,  TJ. OD ROKU -100100
C  KALENDARZA GREGORIANSKIEGO [J1G0=0] ALBO JULIANSKIEGO [J1G0=1].

      J = 4*JD+139361631
      IF(J1G0.LE.0)       J = J+(4*JD+183187720)/146097*3/4*4-3908
      I = MOD(J,1461)/4*5+308
      N = MOD(I,153)/5+1
      M = MOD(I/153,12)+1
      L = J/1461-100100+(8-M)/6
          END

