	Program KH_TeX

c 2010 10 07: poprawka w algorytmie drukowania gorowania Ksiezyca
c      (rzadko, ale pomijal dane); liczne zmiany fornmatu wydruku TeX,
c      gorowanie Cas A zastapione gorowaniem centrum Galaktyki. 
c 2010.08.23-.09.03: Uaktualnienie dtuts do 2014 wg AA 2011
c             kolorowy k_TeX - odroznianie kolorem czasu letniego/zim.
c             dodana tabelka dni swiatecznych na koncu (tylko TeX)
c 2009.10.03: Uaktualnienie dtuts do 2012 wg AA 2009; usuniêcie rysowania
c             dla Surfera.
C 2006 09 13: w glowce dodane minus w h<-18 stop.
	character*1 c
1	write(*,'(a)') ' Kalendarz czy Horoskop [K/k(TeX)/H /q] ? '
	read(*,'(a1)') c
	if(c.eq.'K') call k
	if(c.eq.'k') call k_tex
	if(c.eq.'h'.or.c.eq.'H') call h
	if(c.eq.'q') stop
	go to 1	
	end

	SUBROUTINE INPUT(IY,IM,NM,ISTR,DLG,SZG,NAZWA,krok)
	CHARACTER*1 YES,nazwa*15,text*3
	krok=2
	write(*,1)
1	FORMAT(1H ///,T20,'PROGRAM TEN TWORZY KALENDARZ ASTRONOMICZNY'
     * //' na okres od 1 do 12 miesiecy zadanego roku i miejsca '
     * ,'i umieszcza w zbiorze o    nazwie xxx?0000.TEX lub .TXT, gdzie
     *xxx? pochodzi z nazwy miejsca (? moze byc   znakiem - roku),'
     * ,' a ?0000 to rok. Program '
     * ,'uzywa wlasnych efemeryd cial Ukladu'/' Slonecznego, ktore w '
     * ,'latach ok. 1650 - 2250 maja dokladnosc rzedu 1''',
     * ' (oprocz'/' Plutona [15'']). Aby uzyskac kalendarz nalezy podac'
     * ,' nastepujace dane:'//'  ROK (np. 1987), '
     * /'  PIERWSZY MIESIAC (np. 2 = luty; -5 oznacza 5 kolejnych lat),'
     *,/'  ILOSC MIESIECY (+/-1 do +/-12 [minus '
     * ,'-> kalendarz bez cz. planetarnej]),'/'  NUMER STR'
     * ,'EFY CZASOWEJ (3 = czas moskiewski, 2 = letni polski,',
     * ' 1 = zimowy [CSE],0 = uniwersalny [UT], -4 = atlantycki '
     *,'[Kanada], -5 do -8 = eastern, central,    mountain i pacific'
     *,' [Kanada, USA]),'/'  DLUGOSC (wschodnia) i SZEROKOSC GEOGRAFICZ'
     *,'NA (w stopniach)  oraz'/'  NAZWE MIEJSCA (do 15 znakow; jesli be
     *z spacji, moze byc bez ''uszu'').')

	CALL DATE(nazwa(1:8))
c	 print*,'today:'//nazwa(1:8)//'|'
c dla F77L v.2.2
	if(ichar(nazwa(7:7)).gt.57) nazwa(7:7)=char(ichar(nazwa(7:7))-10)
c dla F77l3 LF90
	read(nazwa,'(i2,1x,i2,1x,i2))') im,nm,iy
	iy=iy+2000
	nm=13-im
	istr=1
	dlg=18.5641
	szg=53.0955
	nazwa='Piwnice, CA UMK'
9	do 10 i=15,2,-1
	if(nazwa(i:i).ne.' ') go to 2
10	continue
2	write(*,3) iy,IM,NM,ISTR,DLG,SZG,NAZWA(1:i)
3	FORMAT(1H /'   CZYTAM DANE; podaj je zgodnie z tym przykladem',
     *' (sam przecinek oznacza powto-'/'  ROK +-M +-nM Str   dlug E   '
     *,'szer N            rzenie wartosci juz wypisanej!)'
     * /1x,i5,1h,,3(i3,1h,),2(f8.4,1h,),1h',a,1h'/)
	read(*,*) IY,IM,NM,ISTR,DLG,SZG,NAZWA
	write(*,6) IY,IM,NM,ISTR,DLG,SZG,NAZWA
6	FORMAT(6X,'Przyjalem nastepujace dane: '/1X,i5,3I3,2F10.5,1X,a15
     * /5X,' O.K.? [ENTER=tak]: ')
	read(*,'(a)') YES
	IF(YES.nE.'Y'.and. YES.nE.'y'.and.yes.ne.' ') go to 9
	if(nm.lt.0) return
	write(*,'(a/a)') ' Czesc planetarna co drugi dzien miesiaca.',
     * '    O.K.? [ENTER=tak, 1 (max. 16) oznacza co 1 (max. 16) dni]:'
	read(*,'(a)') text
	if(text.eq.' '.or.text.eq.'Y'.or.text.eq.'y') return
	read(text,'(i2)') krok
	END

	subroutine k_TeX
	CHARACTER*1 today*8,miejsc(15),ifaz*5,cstr*10,czasy(-12:12)*8
     *,MIES(12)*12,PORY(4)*8,IX(7)*4,fname*12,kolorf(6)*6
	integer*2 mie,idz,irok
	real*8 DJE2,RA,D2
	DIMENSION IW(8),W(4),G(3),F(9)
     * ,WS(9),ZA(9),ALF(9),DEL(9),IFAZ(5)
	DATA PORY/'Wiosna','Lato','Jesie/n','Zima'/, today/' '/
	DATA IX/4h Nie,4h Pon,4h Wto,4h/Sro,4h Czw,4hPi/a,4h Sob/
     * ,IFAZ/4HN/ow ,4h 1kw,5hPe/ln,4h 3kw,4h    /,jf/0/,
     * MIES/8HSTYCZE/N,8HLUTY    ,8HMARZEC  ,8HKWIECIEã,8HMAJ     ,
     * 8HCZERWIEC,8HLIPIEC  ,8HSIERPIEã,8HWRZESIEã,12HPA/XDZIERNIK ,
     * 8HLISTOPAD,8HGRUDZIEã/
     *,czasy/'UT$-$12)','UT$-$11)','UT$-$10)','UT$-$9)','Pacif.)',
     *'Mount.)','Central)','Eastern)','Atlant.)','UT$-$3)','UT$-$2)',
     *'UT-1)','UT)','CSE)','CWE)','Mosk.)','UT+4)','UT+5)',
     *'UT+6)','UT+7)','UT+8)','UT+9)','UT+10)','UT+11)','UT+12)'/

	TRANS(X)=2.2918312*X+.4*FLOAT(IFIX(X*3.819719))
	CALL INPUT(IY,M00,MD,ISTRE,DLD,FID,MIEJSC,krok)

c Formowanie nazwy pliku wyjsciowego
	m=0
	do 201 i=1,4
208	m=m+1
	if(m.gt.15) go to 204
	ic=ichar(miejsc(m))
	if(ic.gt.96.and.ic.lt.123) ic=ic-32	!duze litery
	if(i.eq.1.and.(ic.lt.64.or.ic.gt.91)) go to 208
	if(ic.lt.64.and.ic.gt.57.or.ic.gt.91.or.ic.lt.48) ic=ichar('_')
201	fname(i:i)=char(ic)
204	if(i.lt.5) fname(i:4)='______'
200	if(iy.ge.0) write(fname(5:12),'(i4.4,4h.tex)') iy
	if(iy.lt.0) write(fname(4:12),'(i5.4,4h.tex)') iy
	open(6,file=fname,status='UNKNOWN')

	M0=M00
	 if(M00.le.0) then
	MD=isign(12,MD)
	M0=1
	 endif
c poprawka TDT-UT
	dt=dtuts(iy+.5)/(3600*24)
	c=dt*3.141593*2
	write(cstr,'(1h(,a)') czasy(istre)
	DL=(DLD-dt*360.9864)*3.141593/180.
	FI=FID*3.141593/180.
	CALL DATE(TODAY)
	if(ichar(today(7:7)).gt.57) then
c dla F77L v.2.2
	 read(today,'(i2,1x,i2,2x,i1))') mie,idz,irok
	 irok=irok+(ichar(today(7:7))-58)*10+2000
	else
c dla F77l3 LF90
	 read(today,'(i2,1x,i2,1x,i2))') mie,idz,irok
	 irok=irok+2000
	endif
c nastepna linia dla kompilatora MICROSOFT v. 5.0
c	if(mie*idz.eq.0) CALL getdat(irok,mie,idz)
	write(6,'((a))')   '\documentstyle[colordvi]{article} \prefixing
     *\nofiles\pagestyle{empty}',' ',
     *'% Na laserowej drukowac z opcja /tr1 (dobrac ponizsze offsety)',
     *' ','\hoffset-.65in\voffset-.7in\textwidth9in\textheight6.7in',
     *'\font\SY=yinit scaled\magstep1\def\K{\hbox{\SY m}}\def\Z{\hbox{\S
     *Y 3}}','\begin{document}\special{landscape}\vspace*{-1em}',
     *'\font\Big=cmr10 scaled 5160\font\Bigb=cmbx10 scaled 4300',
     *'\def\RGB#1{\special{color push rgb #1}}',
c     *'\def\blu{\RGB{0. 0. .5}}\def\braz{\RGB{.35 0. 0.}}', ! do wydrukow
     *'\def\blu{\RGB{0. 0. .35}}\def\braz{\RGB{.30 0. 0.}}',
     *'\def\czern{\RGB{.0 0. 0.}}\def\ziel{\RGB{.0 0.15 0.}}'
	WRITE(6,10) iy,miejsc
     *,DLD,FID,ISTRE,cstr,IDZ,MIE,IROK-(irok/100)*100+2000
	CALL SEASON(IY,ISTRE,IW,W)
	MM=NWIELK(IY)
	WRITE(6,112) IY,MIES(MOD(MM,10)),MM/10,nint(dt*(24*3600))
112	FORMAT(1H /1X,'\bigskip\bigskip\noindent Data Wielkanocy: {\bf ',
     * i6,1x,a8,i3,3X,'}\hfill Dla ca/lego roku przyj/eto: TDT $-$ UT ='
     *,i5,2h s/)
	WRITE(6,111) (PORY(I),IW(I),(W(I)-dt*24-amod(W(I)-dt*24,1.)*.4),
     *I=1,4)
111   FORMAT(1H ,'\par\czern\noindent Pory roku [miesi/ac dzie/n,godz.mm
     *]\blu\ (czas zimowy):\hfill{ '/a6,'\bf\ III',I3,',',F6.2,',
     *\rm\hfill ',A4,'\bf\ VI',I3,1H,,F6.2,',\rm\hfill ',A7,'\bf\ IX',
     *I3, 1H,,F6.2,',\rm\hfill ',A4,'\bf\ XII',I3,1h,,F6.2,7h}\czern)

	LetniCzas=1
		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)
c	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)
c	djOct=jdOct-.5d0 + 1d0/24	! 1:00 UTC (= 2:00 CET)

      write(6,'(/a,2(i3,a),i5,a)')'\medskip\braz\noindent Dane w czasie  
     &letnim odr/o/zniono kolorem br/azowym; czas ten obowi/azuje {\bf  
     &od',NMar,' marca do',NOct,' pa/xdziernika',iy,' r.}\czern\\[3mm]'
		endif

	WRITE(6,11) IY
	AJ0=jD(IY,1,0,0)
c	RA=6.113816+(IY-1950)*1.966E-4	! Cassiopea A
c Galactic Center RA 17h45m40.04s, Dec -29° 00' 28.1" (J2000 epoch)
	DJE2=(IY-2000)*365.25d0+2451545.d0
	call PREnew(2451545.d0,DJE2,4.6498509244d0,0.5062817165d0,RA,D2)

	MM=12
	IF(M0.EQ.0) GO TO 55
	M0=IABS(M0)-1
	MM=MOD(MAX0(1,M0+IABS(MD))-1,12)+1
55	M0=M0+1
	F(9)=9.


C  PETLA MIESIECY (MM < 13 (MIESIECY))
	Letni=0
	if(M0.gt.3.and.M0.lt.11) then
	 Letni=-1
	else
	 write(6,'(a)')  '\blu '
	endif
	DO 30 M=M0,MM

300	if(M.ge.3.and.M.lt.11.and.Letni.eq.-1)
     &  write(6,'(a)') '\braz '
     	ISTR=ISTRE+iabs(Letni)

c***	IF(MD.GE.0.OR.MOD(M+1-M0,2).EQ.1) 
	IS=12-ISTR
	WRITE(6,"(1H /,'{\break\noindent\Large\bf\ \,',A12,' \,',i4,
     & '\\[6mm]')") MIES(M),iy

	if(M.ne.3.and.M.ne.10) 	WRITE(6,1) MIEJSC,IS
1	FORMAT('}\break\noindent\halign{\hfil#&\hfil#&\ #\hfil&&',
     *'\hfil#\cr',3('\span')/,15A1,'\hfill&',9('\span')/
     *'\bf\ \leftarrowfill\ Dane na godz.\ 12 czasu strefowego (='
     *,I3,'$^h$ UT) \rightarrowfill\ &\bf\ \dotfill'/
     *6('\span'),'\ Czas strefowy [hhmm.]\ \dotfill\cr'/)
	if(M.eq.3) WRITE(6,41) MIEJSC,IS,IS-1
41	FORMAT('}\noindent\halign{\hfil#&\hfil#&\ #\hfil&&',
     *'\hfil#\cr',3('\span')/,15A1,'\hfill&',9('\span')/
     *'\bf\ \leftarrowfill\ Dane na godz.\ 12 czasu strefowego (=',I3,
     *'\braz\ lub',i3,'$^h$ \blu\ UT) \rightarrowfill\ &\bf\ \dotfill'
     */6('\span'),'\ Czas strefowy [hhmm.]\ \dotfill\cr'/)
	if(M.eq.10) WRITE(6,42) MIEJSC,IS,IS+1
42	FORMAT('}\noindent\halign{\hfil#&\hfil#&\ #\hfil&&',
     *'\hfil#\cr',3('\span')/,15A1,'\hfill&',9('\span')/
     *'\bf\ \leftarrowfill\ Dane na godz.\ 12 czasu strefowego (=',I3,
     *'\blu\ lub',i3,'$^h$ \braz\ UT) \rightarrowfill\ &\bf\ \dotfill'
     */6('\span'),'\ Czas strefowy [hhmm.]\ \dotfill\cr'/)
	AJD=Jd(IY,M,1,0)/1000
	WRITE(6,2)  INT(AJD)! MIES(M),iy,
	MAX=MOD(M/8+M,2)+30
	IF(M.EQ.2)  MAX=28+ (4-MOD(IY,4))/4-
     *        (100-MOD(IY,100))/100+(400-MOD(IY,400))/400
	IFA=0
C  PETLA DNI - CZESC PODSTAWOWA KALENDARZA
	DO 5  ND=1,MAX+1
	IF(MOD(ND,5).EQ.1.AND.MOD(ND,30).NE.1) WRITE(6,8)
	AJ=Jd(IY,M,ND,0)
	ID=AMOD(AJ+1.,7.)+.5
c SUNSET(AJ,DLUGWS,SZEROK,ISTREF,RCZ,WSCH,POLUDN,ZACH,RA,DEKL,RS,DLU)
302	CALL SUNSET(AJ,DL,FI,ISTR,F(2),W(1),GS,W(2),F(3),F(4),RS,DLS)
	CALL MOONSE(AJ,DL,FI,ISTR,F(1),W(3),GK,W(4),F(6),F(7),F(8),DLK,SZ)
	f(1)=f(1)+c
	GOR=DMOD(35.93447d0+(RA-F(1))*3.80929d0,23.93447d0)+8.3333E-3
	IGH=GOR
	IGM=60.*AMOD(GOR,1.)
	F(5)=GS-c
	JR=SIN((DLK-DLS)/2.)**2*100.+.5
	write(ifaz(5),'(i4)') jr
	F9=(DLK-DLS)*.6366198+4.
	F(4)=F(4)*15.
	F(7)=F(7)*15.
	DO 20 I=1,7
	F(I)=TRANS(F(I))
20	IF(I.LT.4.OR.I.EQ.6) F(I)=F(I)-.004*AMOD(F(I)*100.,1.)
	F(2)=F(2)*100.
	JF=AMOD(F9,4.)+1.
	F9=AMOD(F9,1.)*29.53059*.25

	IF(F9.GT.F(9)) JF=5
	  if(nd.le.max) F(9)=F(9)+1.
	IF(JF.EQ.5) GO TO 332

	CALL FAZA(AJ,12.-ISTR,F9,JF)
	DD=ND-F9+.5
	  if(nd.le.max) f(9)=f9
	IF(DD.LT.0.) GO TO 332
	IFA=IFA+1
	ALF(IFA)=IFIX(DD)
		if(ISTR.ne.ISTRE) kolorf(IFA)='\braz '
		if(ISTR.eq.ISTRE) kolorf(IFA)='\blu '
	DD=AMOD(DD-dt,1.)*24.+1./120.
	DEL(IFA)=DD-AMOD(DD,1.)*.4- .005
	WS(IFA)=JF
332	JR=AJ-AJ0
	AJD=AMOD(AJ,1000.)
	G(1)=F(5)*100.
	F(5)=RS
	G(3)=TRANS(GK-c)*100.
	IF(ABS(G(3)-1200.).GT.1200.) G(3)=G(3)+20000.
	DO 21 J=1,4
	WP=TRANS(W(J)+2.181662E-3-c)
	IW(J)=WP
	IW(J+4)=AMOD(WP,1.)*100
	IF(ABS(WP-12.).LE.12.+c*12/3.141593) GO TO 21
	IW(J)=IW(J)+1000
	IW(J+4)=IW(J+4)+1000
21	CONTINUE
	jf=5	! 5 xii 1997 - wywalone tekstowe opisy faz
	g3=nint(G(3)*10)/10.
	gj=mod(G3,10.)
	jg=int(mod(G3,100.))/10
	if(abs(G3).gt.2400.) gj=100
	if(abs(G3).gt.2400.) jg=100
	if(nd.le.max) WRITE(6,100) JR,ND,IX(ID+1),IFIX(AJD),F,IFAZ(JF),
     * (IW(J),IW(J+4),G(J),iw(j+1),IW(J+5),J=1,1),
     * (IW(J),IW(J+4),int(G3/100.),jg,gj,iw(j+1),IW(J+5),
     * J=3,3) ,IGH,IGM ,ND
100	FORMAT(I4,1h&,I3,1h&,A4,1h&,I3.3,1h&,F10.5,1h&,F7.2,1h&,
     * F8.4,1h&,F7.2,1h&,F8.5,1h&,F8.4,1h&,F7.2,1h&,
     * F7.3,1h&,F4.1,1h&,A5,1h&,(I2,i2.2,1h&,F6.1,1h&,I2,i2.2,1h&),
     *(I2,i2.2,1h&,i2,i1,F3.1,1h&,I2,i2.2,1h&)
     *,I2,i2.2,1h&,I3,'\cr')
	if(M.eq.3.and.ND.eq.NMar.and.Letni.ne.-1) then
		Letni=-1
		ISTR=ISTRE+1
		write(6,'(a)') '\braz'
		go to 302
	endif
	if(M.eq.10.and.ND.eq.NOct.and.Letni.ne.0) then
		Letni=0
		ISTR=ISTRE
		write(6,'(a)') '\blu'
		go to 302
	endif
	ID=ID+1
5	CONTINUE
	if(ALF(IFA).gt.MAX) IFA=IFA-1
		if(M.ne.3.and.M.ne.10) then
	WRITE(6,334) (IFAZ(IFIX(WS(J))),IFIX(ALF(J)),DEL(J),J=1,IFA)
334	FORMAT(1H}/'\par\bigskip Fazy Ksi/e/zyca [Dm, h.mm]:',
     &5('\qquad',a8,i3,1h,,f6.2))
		else
	WRITE(6,335) (kolorf(J),IFAZ(IFIX(WS(J))),IFIX(ALF(J)),DEL(J)
     &,J=1,IFA)
335	FORMAT(1H}/'\par\bigskip\czern Fazy Ksi/e/zyca [Dm, h.mm]:',
     &5('\qquad',a6,a8,i3,1h,,f6.2))
	
		endif
	 IF(MD.LT.0) GO TO 30


C  Czesc planetarna kalendarza
	if(M.eq.3) then
		write(6,'(a)') '\blu'
		Letni=0
		ISTR=ISTRE
	endif
	if(M.eq.10) then
		write(6,'(a)') '\braz'
		Letni=-1
		ISTR=ISTRE+1
	endif
	WRITE(6,"(1H /,'\break\noindent{\Large\bf\ \,',A12,' \,',i4,
     & ' }(c.d.)\\[4mm]')") MIES(M),iy
	WRITE(6,12) miejsc 	!mies(m),iy,

	DO 105 ND=1,MAX,krok
	IF(MOD(ND,6).EQ.1) WRITE(6,8)
	AJ=JD(IY,M,ND,0)
122	CALL PLASET(AJ,DL,FI,ISTR,WS,ZA,ALF,DEL,S0)
	S0=S0*12./3.141593
	DO 121 I=1,9
	ALF(I)=ALF(I)*12./3.141593
	DEL(I)=DEL(I)*180./3.141593
	WS(I)= TRANS(WS(I)-c+2.181662E-3)
121	ZA(I)= TRANS(ZA(I)-c+2.181662E-3)
	if((M.eq.10.and.ND.eq.NOct/2*2+1.and.Letni.eq.-1.and.ND.ne.NOct)
     & .or.(M.eq.3.and.ND.eq.NMar/2*2+1.and.Letni.eq.0.and.ND.ne.NMar))
     & go to 124
	WRITE(6,101) ND,(WS(I),ZA(I),I=1,9),(ALF(I),DEL(I),I=1,9),S0,ND
124	if(M.eq.3.and.ND.eq.NMar/2*2+1.and.Letni.eq.0) then
		Letni=-1
		ISTR=ISTRE+1
		write(6,'(a)') '\braz'
		go to 122
	endif
	if(M.eq.10.and.ND.eq.NOct/2*2+1.and.Letni.eq.-1) then
		Letni=0
		ISTR=ISTRE
		write(6,'(a)') '\blu'
		go to 122
	endif
105	continue

	write(6,'(1h},a)') '}\centerline{\box0}'!!!!!'\blu'
	if(M.eq.10) then
		Letni=0
		ISTR=ISTRE
		write(6,'(a)') '\blu'
	endif
	
30	CONTINUE
	call Rocznik(iy)
	write(6,'(/a/)') '\end{document}'
	close(6)
	WRITE(*,'(A)') ' *** O.K. ***  Patrz plik ===> '//fname
	if(M00.ge.-1) return
	M00=M00+1
	iy=iy+1
	go to 200	
2	FORMAT('\span\span\span\cleaders\hbox{$\cdot$}\hfill\ Dzie/n\ 
     *\cleaders\hbox{$\cdot$}\hfill& ST \hfil&$\Delta\alpha$
     *&$\alpha_\odot$\hfil&\ \ $\delta_\odot$\hfil&$R_\odot$\hfil
     *&'/'$\alpha_{\K}$\hfil&$\ \ \delta_{\K}$\hfil&$R_{\K}$\ &',
     *'\span Faza \K\ & \span S/lo/nce\span\hfil&\span Ksi/e/zyc\span &'
     *,'Gal.&',
     *'\cr'/'\ Dr&Dm&Dt &   JD&\,\ h.mmss &\quad\,m.ss&
     *\,h.mmss&\quad\,$^\circ.',1h','$\,\enskip&'/'\quad j astr\ &\ \ 
     *h.mmss&\quad\,$^\circ.',1h','$\,\enskip&\ \quad$\rho_{\Z}$\ &'/
     *'dni&\%\ &Wsch&\ \ G/or\ \ &\ Zach&\ Wsch&\ G/or\ \,&\ Zach&\ G/or
     *&Dm\cr     &   &    &',i6,'\cr')
c4	FORMAT(1H //)
6	FORMAT(1H} /)
c7	FORMAT(1H //////)
8	FORMAT('\noalign{\vskip1ex}')
10    FORMAT(1H /'\centerline{\fbox{\strut\Bigb\Brown{Kalendarz astronom
     *iczny',i5,'}}}\bigskip\bigskip\medskip\vfil\blu'
     */'\centerline{\Big\ \uppercase{',15a1,1x,'}}\bigskip\bigskip\vfil'
     *    /'\czern\centerline{\large\sf Wsp/o/lrz/edne geograficzne:'
     */f10.4,'$^\circ$E',F9.4,'$^\circ$N, Nr strefy czasowej 
     * ($>0$, gdy na wsch/od) =',i3,1x,a8,'}\medskip'
     *// 1X,'\noindent Kalendarz jest gregoria/nskim i jest ',
     *  'oparty na efemerydach Uk/ladu S/lonecznego obliczanych z dok/la
     *dno/sci/a 1'''/       '(S/lo/nce 2", Ksi/e/zyc 10", Pluton 15'').  
     * Czas gwiazdowy ma dok/ladno/s/c oko/lo 0.1 s.  Wschody i zachody' 
     */'obliczono z uwzgl/ednieniem horyzontalnej refrakcji  astronomicz
     *nej (34''), rozci/ag/lo/sci S/lo/nca lub Ksi/e/zyca '/
     *'i paralaksy horyzontalnej Ksi/e/zyca.'/ 'Dla dat odleg/lych o wi/
     *ecej ni/z ok.\ 300 lat od roku 1950 dok/ladno/sci po/lo/ze/n plane
     *t s/a gorsze ni/z 1''!'/'\hfill{\small\it (Kalendarz ten opracowan
     *o [KMB] w ','Katedrze Radioastronomii CA UMK, Toru/n---Piwnice.'/
     *'\quad Oblicz.\ dnia',i3,1h.,i2.2,1h.,i4,' r.)}')
11	FORMAT(1H /' \bigskip\smallskip',
     * '{\it Oznaczenia:}\\[-1.ex]'/'\halign{&#\hfil\cr'/
     * 'Dr &--- & dzie/n',i5,' roku (1 stycznia Dr = 1)&',
     * '$\delta$ &--- & deklinacja (S/lo/nca lub Ksi/e/zyca)\cr'/
     * 'Dm &---& dzie/n miesi/aca&',
     *'$R$ &---& odleg/lo/s/c (S/lo/nca lub Ksi/e/zyca)\cr'/
     * 'Dt &---& dzie/n tygodnia&',
     * 'j astr &---& jednostka astronomiczna\cr'/
     * 'JD &---& data julia/nska modulo 1000 o 12$^h$ UT&',
     * '$\rho_{\Z}$ &---& promie/n Ziemi\cr'/
     * '& &(po nag/l/owku: Entier(JD\normalslash1000) dla Dm = 1)&',
     *'Faza \K&---& dni po kwadrze i \%\ o/swietlonej cz/e/sci tarczy Ks
     *i/e/zyca\cr'/'ST &---& prawdziwy czas gwiazdowy&',
     * 'Wsch &---& wsch/od (S/lo/nca, Ksi/e/zyca lub planety)\cr'/
     * '$\Delta\alpha$&---& r/ownanie czasu (czas s/loneczny prawdziwy',
     * ' minus /sredni)\hskip6em&G/or &---& g/orowanie czyli tranzyt',
     * ' (Gal. = centrum Galaktyki)\cr'/
     * '$\alpha$ &---& rektascensja (S/lo/nca lub Ksi/e/zyca)&',
     * 'Zach &---& zach/od (S/lo/nca, Ksi/e/zyca lub planety)\cr}'//
     * '\vfil\tabskip=.7ex\baselineskip=2.3ex ',
     * '\def\minus{-}\catcode`\-=\active\def-{$\minus$}')
12	FORMAT(1H /
     *'\begin{center}\large\sf/Swity, zmierzchy astronomiczne ($h_\odot=  
     * \minus18^\circ$; kolumna S/LO/NCE), wschody i zachody planet  
     * (w czasie strefowym [hh.mm])\\'/'oraz, w drugim wierszu dla danej 
     * daty, wsp/o/lrz/edne (rektascensja [godziny] i deklinacja [stopni 
     *e]) S/lo/nca i planet\\'/'w po/ludnie strefowe dla miejsca ',15a1,
     *' oraz czas gwiazdowy o p/o/lnocy UT w Greenwich [godziny]'/
     *'\rm\end{center}\medskip'//
     *'\setbox0=\vbox{\noindent\halign{&\hfil#\cr'/
     *' &\span\bf S/LO/NCE \hfil&\span\bf MERKURY&\span\bf WENUS 
     *\hfil&\span\bf MARS \hfil&\span\bf JOWISZ'/'\hfil&\span\bf SATURN
     *\hfil&\span\bf URAN \hfil&\span\bf NEPTUN \hfil&\span\bf PLUTON
     *\hfil&\bf ST \hfil&\cr'/'Dm& /Swit & Zm',3h'ch,8('& Wsch & Zach')
     *,'&(0$^h$UT)&Dm\cr'/'\noalign{\vskip.7ex}')
101	FORMAT(I3,1h&,9(F5.2,1H&,F5.2,1h&),'\cr'/
     *1X,9(1h&,F6.3,1h&,F6.2),1h&, F9.5,1h&,I3,'\cr')
	END

	subroutine K

C                     KALENDARZ ASTRONOMICZNY

C  PROGRAM OBLICZA EFEMERYDY UKLADU SLONECZNEGO Z DOKLADNOSCIA OKOLO
C  1 MINUTY LUKU (15' DLA PLUTONA) I CZAS GWIAZDOWY (PRAWDZIWY) Z DOK-
C  LADNOSCIA OK. 0.1 SEK.    WSCHODY I ZACHODY OBLICZANE SA Z UWZGLED-
C  NIENIEM REFRAKCJI ASTRONOMICZNEJ (34'), ROZCIAGLOSCI SLONCA LUB
C  KSIEZYCA I PARALAKSY HORYZONTALNEJ KSIEZYCA.    JAKO DANE POCZATKOWE
C  NALEZY PODAC (W JEDNYM REKORDZIE NA KAZDY KALENDARZ): ROK (SPOMIEDZY
C  OK. R. 1650 I 2250), MIESIAC PIERWSZY (1 DO 12 (0, JESLI CALY ROK)),
C  ILOSC MIESIECY (1 DO 12 LUB -1 DO -12 (UJEMNA WARTOSC OZNACZA KALEN-
C  DARZ BEZ CZESCI PLANETARNEJ)), NUMER STREFY CZASOWEJ (W PRAKTYCE OD
C  -12 (NA ZACHODZIE) DO +13 (NA WSCHODZIE)),
C  DLUGOSC (WSCHODNIA) I SZEROKOSC GEOGRAFICZNA (W STOPNIACH)
C  ORAZ NAZWE MIEJSCOWOSCI LUB MIEJSCA (CO NAJWYZEJ 15 ZNAKOW)

C   WYKORZYSTYWANE NIESTANDARTOWE PODPROGRAMY I FUNKCJE:

C   FPAGE, INPUT, SUNSET, MOONSE, PLASET, PLANET, FAZA, SUN,
C   KSIEZY, MERKUR, WENUS, MARS, JOWISZ, SATURN, URAN, NEPTUN, PLUTON,
C   SEASON, JD, NWIELK, DATE

	CHARACTER*1 today*8,miejsc(15),kl(9),na*9,cstr*10,fname*12
     *,czasy(-12:12)*8,      MIES(12)*8,PORY(4)*8,IX(7)*4,IFAZ(5)*5
	integer*2 mie,idz,irok
	DIMENSION IW(8),W(4),G(3),F(9)
     * ,WS(9),ZA(9),ALF(9),DEL(9)
	DATA PORY/'WIOSNA','LATO','JESIEN','ZIMA'/, today/'        '/
	DATA IX/4H NIE,4H PON,4H WTO,4H SRO,4H CZW,4H PIA,4H SOB/
     * ,IFAZ/4H NOW ,4H 1KW,4HPELN,4H 3KW,4H    /,JF/0/,
     * MIES/8HSTYCZEN ,8HLUTY    ,8HMARZEC  ,8HKWIECIEN,8HMAJ     ,
     * 8HCZERWIEC,8HLIPIEC  ,8HSIERPIEN,8HWRZESIEN,8HPAZDZ.  ,8HLISTOPAD
     * ,8HGRUDZIEN/,KL/1HK,1HA,1HL,1HE,1HN,1HD,1HA,1HR,1HZ/
     *,czasy/'(UT-12)','(UT-11)','(UT-10)','(UT-9)','(Pacif.)',
     *'(Mount.)','(Centr.)','(East.)','(Atlant)','(UT-3)','(UT-2)',
     *'(UT-1)','(UT)','(CSE)','(CWE)','(Mosk.)','(UT+4)','(UT+5)',
     *'(UT+6)','(UT+7)','(UT+8)','(UT+9)','(UT+10)','(UT+11)','(UT+12)'/

	TRANS(X)=2.2918312*X+.4*FLOAT(IFIX(X*3.819719))
	CALL INPUT(IY,M00,MD,ISTR,DLD,FID,MIEJSC,krok)

c Formowanie nazwy pliku wyjsciowego
	m=0
	do 201 i=1,4
208	m=m+1
	if(m.gt.15) go to 204
	ic=ichar(miejsc(m))
	if(ic.gt.96.and.ic.lt.123) ic=ic-32	!duze litery
	if(i.eq.1.and.(ic.lt.64.or.ic.gt.91)) go to 208
	if(ic.lt.64.and.ic.gt.57.or.ic.gt.91.or.ic.lt.48) ic=ichar('_')
201	fname(i:i)=char(ic)
204	if(i.lt.5) fname(i:4)='______'
200	if(iy.ge.0) write(fname(5:12),'(i4.4,4h.txt)') iy
	if(iy.lt.0) write(fname(4:12),'(i5.4,4h.txt)') iy
	open(6,file=fname,status='UNKNOWN')

	M0=M00
	 if(M00.le.0) then
	MD=isign(12,MD)
	M0=1
	 endif
c poprawka TDT-UT
	dt=dtuts(iy+.5)/(3600*24)
	c=dt*3.141593*2
	cstr=czasy(istr)
	WRITE(6,4)
58	WRITE(6,7)
	CALL FPAGE(KL,9)
	WRITE(6,7)
	write(na,'(2hNA,i5.4,2h R)') iy
	CALL FPAGE(NA,9)
	WRITE(6,7)
	CALL FPAGE(MIEJSC,9)
	DL=(DLD-dt*360.9864)*3.141593/180.
	FI=FID*3.141593/180.
	CALL DATE(TODAY)
	if(ichar(today(7:7)).gt.57) then
c dla F77L v.2.2
	 read(today,'(i2,1x,i2,2x,i1))') mie,idz,irok
	 irok=irok+(ichar(today(7:7))-58)*10+2000
	else
c dla F77l3 LF90
	 read(today,'(i2,1x,i2,1x,i2))') mie,idz,irok
	 irok=irok+2000
	endif
c nastepna linia dla kompilatora MICROSOFT v. 5.0
c	if(mie*idz.eq.0) CALL getdat(irok,mie,idz)
	WRITE(6,10) DLD,FID,ISTR,cstr,IDZ,MIE,IROK-(irok/100)*100+2000
	CALL SEASON(IY,ISTR,IW,W)
111	FORMAT(1H ,4X,'PORY ROKU [MIES. DZIEN,GODZ.MM]:  ',A6,' III',I3,
     * ',',F6.2,'   *** ',A4,' VI',I3,1H,,F6.2,'   *** ',A6,' IX',I3,
     * 1H,,F6.2,'   *** ',A4,' XII',I3,1h,,F6.2)
	MM=NWIELK(IY)
	WRITE(6,112) IY,MIES(MOD(MM,10)),MM/10,nint(dt*(24*3600))
112	FORMAT(1H /5X,'DATA WIELKANOCY: ',I6,1X,A8,I3,50x,
     * 'Przyjeto dla calego roku: TDT - UT =',i5,2h s)
	WRITE(6,111) (PORY(I),IW(I),(W(I)-dt*24-amod(W(I)-dt*24,1.)*.4),
     *I=1,4)
	WRITE(6,11) IY
	AJ0=jD(IY,1,0,0)
	RA=6.113816+(IY-1950)*1.966E-4
	MM=12
	IF(M0.EQ.0) GO TO 55
	M0=IABS(M0)-1
	MM=MOD(MAX0(1,M0+IABS(MD))-1,12)+1
55	M0=M0+1
	IS=12-ISTR
	F(9)=9.
C  PETLA MIESIECY (MM < 13 (MIESIECY))
	DO 30 M=M0,MM
	IF(MD.GE.0.OR.MOD(M+1-M0,2).EQ.1) WRITE(6,1) MIEJSC,IS
	AJD=Jd(IY,M,1,0)/1000
	WRITE(6,2) IY,MIES(M)
	WRITE(6,60) INT(AJD)
	MAX=MOD(M/8+M,2)+30
	IF(M.EQ.2)  MAX=28+ (4-MOD(IY,4))/4-
     *        (100-MOD(IY,100))/100+(400-MOD(IY,400))/400
	IFA=0
C  PETLA DNI - CZESC PODSTAWOWA KALENDARZA
	DO 5  ND=1,MAX+1
	IF(MOD(ND,5).EQ.1.AND.MOD(ND,30).NE.1) WRITE(6,8)
	AJ=Jd(IY,M,ND,0)
	ID=AMOD(AJ+1.,7.)+.5
c SUNSET(AJ,DLUGWS,SZEROK,ISTREF,RCZ,WSCH,POLUDN,ZACH,RA,DEKL,RS,DLU)
	CALL SUNSET(AJ,DL,FI,ISTR,F(2),W(1),GS,W(2),F(3),F(4),RS,DLS)
	CALL MOONSE(AJ,DL,FI,ISTR,F(1),W(3),GK,W(4),F(6),F(7),F(8),DLK,SZ)
c            MOONSE(AJ,DLG,SZG,ISTR,SID,WSCH,TK,ZACH,RAK,DEK,RK,DLK,SZK)
	f(1)=f(1)+c
	GOR=AMOD(12.+(RA-F(1))*3.80929,23.93447)+8.3333E-3
	IGH=GOR
	IGM=60.*AMOD(GOR,1.)
	F(5)=GS-c
	JR=SIN((DLK-DLS)/2.)**2*100.+.5
	write(ifaz(5),'(i4)') jr
	F9=(DLK-DLS)*.6366198+4.
	F(4)=F(4)*15.
	F(7)=F(7)*15.
	DO 20 I=1,7
	F(I)=TRANS(F(I))
20	IF(I.LT.4.OR.I.EQ.6) F(I)=F(I)-.004*AMOD(F(I)*100.,1.)
	F(2)=F(2)*100.
	JF=AMOD(F9,4.)+1.
	F9=AMOD(F9,1.)*29.53059*.25

	IF(F9.GT.F(9)) JF=5
	  if(nd.le.max) F(9)=F(9)+1.
	IF(JF.EQ.5) GO TO 332

	CALL FAZA(AJ,12.-ISTR,F9,JF)
	DD=ND-F9+.5
	  if(nd.le.max) f(9)=f9
	IF(DD.LT.0.) GO TO 332
	IFA=IFA+1
	ALF(IFA)=IFIX(DD)
	DD=AMOD(DD-dt,1.)*24.+1./120.
	DEL(IFA)=DD-AMOD(DD,1.)*.4- .005
	WS(IFA)=JF
332	JR=AJ-AJ0
	AJD=AMOD(AJ,1000.)
	G(1)=F(5)*100.
	F(5)=RS
	G(3)=TRANS(GK-c)*100.
	IF(ABS(G(3)-1200.).GT.1200.) G(3)=G(3)+20000.
	DO 21 J=1,4
	WP=TRANS(W(J)+2.181662E-3-c)
	IW(J)=WP
	IW(J+4)=AMOD(WP,1.)*100
	IF(ABS(WP-12.).LE.12.+c*12/3.141593) GO TO 21
	IW(J)=IW(J)+1000
	IW(J+4)=IW(J+4)+1000
21	CONTINUE
	ID=ID+1
	jf=5	! 5 xii 1997 - wywalone tektowe opisy faz
	if(nd.le.max) WRITE(6,100) JR,ND,IX(ID),IFIX(AJD),F,IFAZ(JF),
     * (IW(J),IW(J+4),G(J),iw(j+1),IW(J+5),J=1,3,2) ,IGH,IGM ,ND
5	CONTINUE
	WRITE(6,334) (IFAZ(IFIX(WS(J))),IFIX(ALF(J)),DEL(J),J=1,IFA)
334	FORMAT(1H /'  FAZY KSIEZYCA [DM,H.MM]:',5(2X,A4,I3,F6.2,:'  ***'))
	 IF(MD.LT.0) WRITE(6,6)
	 IF(MD.LT.0) GO TO 30
C  CZESC PLANETARNA KALENDARZA
	WRITE(6,12)
	DO 105 ND=1,MAX,krok
	IF(MOD(ND,6).EQ.1) WRITE(6,8)
	AJ=JD(IY,M,ND,0)
	CALL PLASET(AJ,DL,FI,ISTR,WS,ZA,ALF,DEL,S0)
	S0=S0*12./3.141593
	DO 121 I=1,9
	ALF(I)=ALF(I)*12./3.141593
	DEL(I)=DEL(I)*180./3.141593
	WS(I)= TRANS(WS(I)-c+2.181662E-3)
121	ZA(I)= TRANS(ZA(I)-c+2.181662E-3)
105	WRITE(6,101) ND,(WS(I),ZA(I),I=1,9),(ALF(I),DEL(I),I=1,9),S0,ND
30	CONTINUE
	close(6)
	WRITE(*,'(A)') ' *** O.K. ***         Patrz plik ===> '//fname
	if(M00.ge.-1) return
	M00=M00+1
	iy=iy+1
	go to 200	
1	FORMAT(1H //,2X,15A1,2X,11(1H*),' DANE NA GODZ. 12 CZASU STREF. (=
     * ',I3,' GODZ. UT) ',10(1H*),3X,'======== CZAS STREFOWY [HHMM.] =='
     * ,'=======',/)
2	FORMAT(1H ,I5,A9,6X,'ST',7X,'R CZ  RA SL  DEK SL R SL     RA KS',
     *  '  DEK KS R KS   FAZA KS        SLONCE            KSIEZYC'
     * ,4X,' CAS A')
4	FORMAT(1H //)
6	FORMAT(1H /)
7	FORMAT(1H //////)
8	FORMAT(1H )
10	FORMAT(1H ///14X,'WSPOLRZEDNE GEOGRAFICZNE:',F10.4,1HE,F9.4,
     * 'N [stopni], NR STREFY CZASOWEJ (>0, gdy na wschod) =',i3,1x,a8//
     * 7X,'Ponizszy kalendarz jest gregorianskim i jest oparty na ',
     * 'efemerydach Ukladu Slonecznego obliczanych z dokladnoscia ',
     *'1 minuty luku',/,5X,'(Slonce 2", Ksiezyc 10", Pluton 15'').   Cza
     *s gwiazdowy ma dokladnosc okolo 0.1 sek.    Wschody i zachody obli
     *czono z uwzglednie-'/5x,'niem horyzontalnej refrakcji astronomiczn
     *ej (34''), rozciaglosci Slonca lub Ksiezyca ',
     *'i paralaksy horyzontalnej Ksiezyca.'/6x,' Dla dat odleglych o wie
     *cej niz ok. 300 lat od roku 1950 dokladnosci polozen planet',
     *' sa gorsze niz 1''!'/5x,'(Kalendarz ten wykonano w '
     * ,'Katedrze Radioastronomii UMK, Torun - Piwnice.   Oblicz. dnia'
     *,i3,1h.,i2.2,1h.,i4,' r.)',20x,'[KB]')
11	FORMAT(1H ,//,5X,
     * 'OZNACZENIA',//,5X,'DR - dzien',i5,' roku (1 stycznia DR=1)',
     * T70,'DEK SL(KS) - deklinacja Slonca (Ksiezyca)'/
     * 5X,'DM - dzien miesiaca',
     * T70,'R SL(KS)   - odleglosc Slonca (Ksiezyca)'/
     * 5X,'DT - dzien tygodnia',
     * T71,'   J ASTR - jednostka astronomiczna'/
     * 5X,'JD - data julianska modulo 1000 w poludnie uniwersalne',
     * T71,'   R ZIEM - promien Ziemi'/
     * 10x,'(po naglowku podano Entier(JD/1000) dla DM=1)',
     *T70,'FAZ KS - dni po kwadrze i procent oswietl. tarczy Ksiezyca'/
     * 5X,'ST - prawdziwy czas gwiazdowy',
     * T70,'WSCH   - wschod (Slonca, Ksiezyca lub planety)'/
     * 5X,'R CZ - rownanie czasu (sloneczny prawdziwy minus sredni)',
     * T70,'GOROW  - gorowanie (przejscie przez poludnik, tranzyt)'/
     * 5X,'RA SL(KS) - rektascensja Slonca (Ksiezyca)',
     * T70,'ZACH   - zachod (Slonca, Ksiezyca lub Planety)'///)
12	FORMAT(1H //4X,'Swity, zmierzchy astronomiczne (h<-18 stop.; kolum
     *na SLONCE), wschody i zachody planet (w czasie strefowym [HH.MM])
     * oraz,'/4x,'w drugim wierszu po dacie, wspolrzedne (rektascensja [
     *godziny] i deklinacja [stopnie]) Slonca i planet w poludnie strefo
     *we'/t40,'oraz czas gwiazdowy o polnocy UT w Greenwich [godziny]',
     *//9X,'SLONCE       MERKURY       WENUS        MARS        JOWISZ'
     * ,T75,'SATURN        URAN        NEPTUN       PLUTON'/
     * '  DM',2X,' SWIT==ZM',3H'CH,8('  WSCH===ZACH'),'  ST(0UT) DM')
101	FORMAT(1H ,I3,1X,9(1X,F5.2,2H==,F5.2),/,5X,9(1X,F6.3,F6.2),
     * F9.5,I3)
60	FORMAT(1H ,'  DR DM DT    JD   H.MMSS    M.SS  H.MMSS  ST.MM J AST
     *R   H.MMSS  ST.MM R ZIEM DNI   %   WSCH  GOROW ZACH  WSCH  GOROW '
     * ,'ZACH GOROW',3H DM/I17)
100	FORMAT(1H ,I4,I3,A4,2X,I3.3,F10.5,F7.2,F8.4,F7.2,F8.5,F8.4,F7.2,
     * F7.3,F4.1,A4,1X,2(2X,2I2,1X,F6.1,1X,2I2),2X,2I2,I3)
	END

	function dtuts(y)
	dimension idt(32),dt(235)		! dt(235=2014-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 - 2014 (uaktualniono 2010.08.19 wg. Astron.Almanac 2011)
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 data extrapolated to 2010 - 2014)
cc    2002   3     4     5     6    7     8     9    10   
     *64.3,64.47,65.57,64.69,64.85,65.15,65.46,65.78,66.1,
cc    2011 12  13 14
     *66.4,67.,67,67./

c Uaktualnienie tabeli dt: dopisanie N nowych wartosci powiazac ze
c zwiekszeniem wymiaru tabeli dt o N oraz w nastepnym wierszu
c powiekszyc graniczny rok po "if(iy.gt." takze o N.

	iytab=2014	!!! -> 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 wartosc skoku 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

	SUBROUTINE PREnew(DJE1,DJE2,ra1,d1,ra2,d2)

C CALCULATES GENERAL PRECESSION FROM DJE1 TO DJE2 using new IAU theory
c Stale wg J.H.Lieske 1979, Astron. Astrophys. 73, 282.
	IMPLICIT REAL*8 (D,r)
	DATA DCSAR/4.848136812D-6/,J2000/2451545/,DC1/2306.2181D0/,
     * DC2/1.39656D0/,DC3/.30188D0/,DC4/.017998D0/,DC5/1.09468D0/
     * ,DC6/2004.3109D0/,DC7/-.8533D0/,DC8/-.42665D0/,DC9/-.041833D0/
	DT0=(Dje1-j2000)/36525
	DT=(Dje2-Dje1)/36525
	DTS=DT*DT
	DTC=DTS*DT
	DZETA=((DC1+(DC2-1.39d-4*dt0)*DT0)*DT+(DC3-3.44d-4*dt0)*DTS+
     * DC4*DTC)*DCSAR
	dzet =DZETA+((DC5-dc3+4.1d-4*dt0)*DTS+.000205*dtc)*DCSAR
	dth=((DC6+(DC7-2.17d-4*dt0)*DT0)*DT+(DC8-2.17d-4*dt0)*DTS+
     * DC9*DTC)*DCSAR

c convert Dj1 position to Dj2 using new IAU1976 precession
	dc=cos(ra1+dzeta)
	ra2=dzet+datan2(sin(ra1+dzeta),-tan(d1)*sin(dth)+cos(dth)*dc)
	d2=asin(sin(d1)*cos(dth)+cos(d1)*sin(dth)*dc)
	ra2=dmod(ra2+6.283185307178d0,6.283185307178d0)
	END

	SUBROUTINE SUNSET(AJ,DLUGWS,SZEROK,ISTREF,RCZ,
     * WSCH,POLUDN,ZACH,RA,DEKL,RS,DLU)
	common /suncom/ xi0,eta0,s,ce,dla
	DOUBLE PRECISION DJD
	DJD=AJ-DBLE(ISTREF/24.)
	CALL Sun(DJD,RA,DEKL,DLU,RS,RCZ,st)
	POLUDN=6.283185*(.5+ISTREF/24.)-DLUGWS-RCZ
	R=(-.0145431-SIN(DEKL)*SIN(SZEROK))/(COS(DEKL)*COS(SZEROK))
	WSCH=1000.
	ZACH=1000.
	IF(ABS(R).GE.1.) RETURN
	R=ACOS(R)
	C=TAN(SZEROK)*COS(RA)*S*.00273791/(SIN(R)*COS(DEKL)**2)
	WSCH=C*3.141593+(POLUDN-R)*(1-C)
	ZACH=-C*3.141593+(POLUDN+R)*(1+C)
	IF(ZACH.GT.6.283185) ZACH =ZACH-6.283185*(1.-C)
	IF(WSCH.LT.0.) WSCH =WSCH+6.283185*(1.+C)
	END

	SUBROUTINE MOONSE(AJ,DLG,SZG,ISTR,SID,WSCH,TK,ZACH,
     *      RAK,DEK,RK,DLK,SZK)
	 COMMON /suncom/xi0,eta0,se,ce,DLA
	DOUBLE PRECISION DJD
	DATA PI,E/3.1415927,.02/
	HA(D)=ACOS(AMIN1(1.,AMAX1(-1.,(H0-SIN(D)*SIN(SZG))/
     * (COS(D)*COS(SZG)))))*.1591549*(2-I)
	UT=.5
	DJD=AJ-2451545D0-ISTR/24.-.5D0
	DO 2 I=1,3,2
	J=-1
1	DJD=DJD+UT
	CALL KSIEZY(DJD+2451545d0,RA,D,DL,B,R)
	ST=3.14159265359D0*DMOD(1.55811454652D0+DMOD(DJD+DJD,2D0)+DJD*
     * (.547581870159D-2+DJD*(1.61549D-15-DJD*1.473D-24))+DLG/PI,2D0)
	DJD=DJD-UT
	H0=SIN(.727507/R-.989E-2)
	DH=1.+.0365*(.917-SIN(D)*SIN(B))/COS(D)**2
	IF(I*J.NE.-1) GO TO 4
	IF(UT.NE..5) GO TO 5
	SID=AMOD(ST+PI+PI+DLA*ce,PI+PI)
	RAK=RA
	DEK=D
	DLK=DL
	SZK=B
	RK=R
5	TK=AMOD((RA-ST)/6.283185+1.,1.)
	IF(TK.GT..5) TK=TK-1.
	IF(ABS(TK).LT.E) GO TO 6
	UT=UT+TK*DH
	GO TO 1
6	TK=(UT+TK*DH)*(PI+PI)
4	H=HA(D)
	DW=AMOD((RA-ST)/6.283185+H+2.,1.)
	IF(DW.GT..5) DW=DW-1.
	D=D+.0917*COS(RA)*DW*DH
	H=HA(D)
	DW=AMOD((RA-ST)/6.283185+H+2.,1.)
	IF(DW.GT..5) DW=DW-1.
	UT=UT+DW*DH
	J=J+ISIGN(1,J)
	IF(ABS(DW).GT.E.AND.IABS(J).LE.3) GO TO 1
	WSCH=UT*(PI+PI)
	IF(ABS(WSCH-PI).LE.PI.OR.J.GT.0) GO TO 3
	J=1
	UT=UT-SIGN(DH,WSCH)
	GO TO 1
3	IF(I.EQ.1) ZACH=WSCH
2	UT=AMOD(UT+(1.-2*H)*DH,DH)

c	print*,' Wsch, zach Ksiez. [h]:', WSCh/pi*12,Zach/pi*12

	IF(ABS(WSCH-ZACH).GT..0005) RETURN
	WSCH=1000.+WSCH
	ZACH=1000.+ZACH
	  END

	SUBROUTINE PLASET(AJ,DL,FI,ISTR,RI,SE,RA,D,S0)
	common /suncom/ xi0,eta0,s,c,dla
	DOUBLE PRECISION DJD,T
	DIMENSION RI(9),SE(9),RA(9),D(9)
	DATA PI/3.141593/,S18,S34/.309017,9.89004E-3/,
     * RK/2.730434E-3/	!,C18/.951057/
	CHA(DD,SH0)=-(SH0+SIN(DD)*SFI)/(COS(DD)*CFI)
	SI(T,DLG)=3.14159265359D0*DMOD(1.55811454652D0+DMOD(T+T,2D0)+T*
     * (.547581870159D-2+T*(1.61549D-15-T*1.473D-24))+DLG/PI,2D0)
	SFI=SIN(FI)
	CFI=COS(FI)
	DJD=AJ-2451545D0-ISTR/24D0
	CALL SUN(DJD+2451545d0,RA(1),D(1),DLU,R,RCZ,st)
	S0=AMOD(SI(AJ-2451545.5D0,0.)+PI+PI+DLA*C,PI+PI)
	ST=AMOD(SI(DJD,DL)+PI+PI+DLA*C,PI+PI)
	DO 1 I=1,9
	RI(I)=1000.
1	SE(I)=1000.
	cH=CHA(D(1),S18)
	IF(ABS(cH).GE.1.) GO TO 2
	  H=ACOS(cH)
CC	  CC=TAN(ASIN((SFI+SIN(D(1))*S18)/(COS(D(1))*C18)))*S*
CC     *   COS(RA(1))*RK/COS(D(1))

c	C=TAN(SZEROK)*COS(RA)*S*.00273791d0/(SIN(R)*COS(DEKL)**2)
c pe³ne wyra¿enie - wazne dla duzych h:
	CC=Rk*COS(RA(1))*S*(sfi/cfi-cH*tan(d(1)))/SIN(H)
c	 =RK*COS(RA)*S*(sin(SZEROK)-sinh*sin(dekl))
c     & /(cos(szerok)*SIN(H)*COS(DEKL)**2)

	  P=PI*ISTR/12.-DL-RCZ
	  RI(1)=PI-(1-CC)*(H-P)
	  SE(1)=PI+(1+CC)*(H+P)
	IF(RI(1).LT.0.) RI(1)=RI(1)+(PI+PI)*(1+CC)
	IF(SE(1).GT.PI+PI) SE(1)=SE(1)-(PI+PI)*(1-CC)
2	CONTINUE
	DO 3 I=2,9
	CALL PLANET(I-1,RA(I),D(I),DLE,SZE,R)
	H=CHA(D(I),S34)
	IF(ABS(H).GT.1.) GO TO 3
	  H=ACOS(H)
	  RI(I)=(AMOD(RA(I)-ST-H+3.*PI,PI+PI)*365.2422+PI)*RK
	  SE(I)=(AMOD(RA(I)-ST+H+3.*PI,PI+PI)*365.2422+PI)*RK
3	CONTINUE
	END

	SUBROUTINE SEASON(IY,ISTR,IDZIEN,GODZ)

C  PROGRAM OBLICZA PRZYBLIZONE (Z DOKL. OK. 1 min. TDT) POCZATKI POR
C  ROKU (IY) W CZASIE STREFOWYM (ISTR JEST NUMEREM STREFY) W DNIACH
C  (W TABLICY IDIEN(I)) I GODZINACH (W GODZ(I), I=1,4 ODPOWIADAJACE
C  WIOSNIE, LATU, JESIENI I ZIMIE ASTRONOMICZNEJ).

	DIMENSION GODZ(4),IDZIEN(4)
	DOUBLE PRECISION DJ
	DO 2 I=1,4
	DD=0.
	DJ=Jd(IY,I*3,20,0)-.5D0-ISTR/24.d0
	DO 1 N=1,6
	CALL Sun(DJ+DD,RA,D,DL,R,RCZ,st)
	DF=I-1-DL/3.14159265d0*2.
	if(abs(df).gt.2.) df=df-sign(4.,df)
	dd=dd+df*91.3105
	if(abs(df).lt.1e-7) go to 3
1	continue
3	IDZIEN(I)=20+DD-DF*45.5
2	GODZ(I)=AMOD(DD-DF*45.5+20.,1.)*24.
	END

	FUNCTION NWIELK(IR)

C  OBLICZA DATE WIELKANOCY ROKU IR.     WYNIK JEST POSTACI N*10+M, 
C  GDZIE N JEST DNIEM MIESIACA M [M = 3 (MARZEC) LUB 4 (KWIECIEN)]

	IMPLICIT INTEGER*4 (A-Z)
	A=MOD(IR,19)
	B=IR/100
	C=MOD(IR,100)
	H=MOD(19*A+B-B/4-(B+1-(B+8)/25)/3+15,30)
	L=MOD(32+2*(MOD(B,4)+C/4)-H-MOD(C,4),7)
	M=(A+11*H+22*L)/451
	N=H+L-7*M+114
	NWIELK=14+10*MOD(N,31)
	IF(N/31.EQ.3) NWIELK=NWIELK-1
	END

	SUBROUTINE FAZA(AJ,UT,DD,KWADRA)
	DOUBLE PRECISION DJ
	DJ=AJ+(UT-12.)/24d0
	DD=0.
1	CALL Sun(DJ-DD,RA,DEC,DS,R,EQT,st)
	CALL KSIEZY(DJ-DD,RA,DEC,DK,SK,R)
	FAZ=AMOD((DK-DS)*2./3.141593+4.,4.)
	IF(KWADRA.LE.0) KWADRA=FAZ+1.
	FAZ=KWADRA-1-FAZ
	IF(DD.EQ.0..AND.FAZ.GT.0.) FAZ=FAZ-4.
	IF(DD.NE.0..AND.ABS(FAZ).GT.2.) FAZ=FAZ-SIGN(4.,FAZ)
	DD=DD-FAZ*29.53059/4.
	IF(ABS(FAZ).GT..00001) GO TO 1
C  O.01 FAZ'A ODPOWIADA OK. 100 MIN. DD
	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

	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

	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 PLANET(NR,RA,DEC,DLE,SZE,R)
	GO TO (1,2,3,4,5,6,7,8), NR
	WRITE (*,'(A,I3)') ' *** Niewlasciwy numer planety:',NR
	RETURN
1	CALL MERKUR(RA,DEC,DLE,SZE,R)
	                            RETURN
2	CALL  WENUS(RA,DEC,DLE,SZE,R)
	                            RETURN
3	CALL   MARS(RA,DEC,DLE,SZE,R)
	                            RETURN
4	CALL JOWISZ(RA,DEC,DLE,SZE,R)
	                            RETURN
5	CALL SATURN(RA,DEC,DLE,SZE,R)
	                            RETURN
6	CALL   URAN(RA,DEC,DLE,SZE,R)
	                            RETURN
7	CALL NEPTUN(RA,DEC,DLE,SZE,R)
	                            RETURN
8	CALL PLUTON(RA,DEC,DLE,SZE,R)
	END

	SUBROUTINE HELGEO(DLH,SZH,RH,RA,DEC,DLE,SZE,R)

C     PODPROGRAM TEN PRZELICZA HELIOCENTRYCZNE WSPOLRZEDNE PLANET (DLH,
C    SZH I RH) NA WSPOLRZEDNE GEOCENTRYCZNE (RA, DEC, DLE, SZE I R)
C  
	COMMON /SUNCOM/ XI0,ETA0,SE,CE,DL
	C = RH*COS(SZH)
	X = C*COS(DLH + DL) + XI0
	C = C*SIN(DLH + DL) + ETA0
	B = RH*SIN(SZH)
	Y = C*CE - B*SE
	Z = C*SE + B*CE
	R   = SQRT(X*X + Y*Y + Z*Z)
	RA  = AMOD(ATAN2(Y,X) + 6.2831853,6.2831853)
	DEC = ASIN(Z/R)
	DLE = AMOD(ATAN2(Y*CE + Z*SE,X) + 6.2831853,6.2831853)
	SZE = ASIN((Z*CE - Y*SE)/R)
	END

	SUBROUTINE MERKUR(RA,DEC,DLE,SZE,R)                            

C    PROCEDURA OBLICZA WSPOLRZEDNE ROWNIKOWE (RA I DEC) I EKLIPTYCZNE
C   (DLE I SZE) [RAD] ORAZ ODLEGLOSC OD ZIEMI [JEDN. ASTR.] MERKUREGO
C   NA MOMENT ODLEGLY O T DNI OD POLUDNIA DNIA 1 STYCZNIA  2000 ROKU. 
C   WYMAGA ONA UPRZEDNIEGO WYWOLANIA PODPROGRAMU SLONCE.

	COMMON /PLAN/ T
c,A1,A7,A8,A13,A16,A19
	DOUBLE PRECISION T,D1,D2
	F(D1,D2) = 6.28318530718D0*DMOD(D1 + T*D2,1D0)
	A9 = F(.700695D0,.011367714D0)
	A10 = F(.485541D0,.01136759566D0)
	A11 = F(.566441D0,.01136762384D0)
	A13 = F(.140023D0,4.45036173D-3)
	DLH = A9+((84386.+2.2E-4*SNGL(T))*SIN(A10)+10733.*SIN(A10+
     * A10)+1892.*SIN(3.*A10)-646.*SIN(A11+A11)+381.*SIN(4.*A10)-306.*
     * SIN(A10-A11-A11)-274.*SIN(A10+A11+A11)-92.*SIN(2.*(A10+A11))+83.*
     * SIN(5.*A10)-28.*SIN(3.*A10+2.*A11)+25.*SIN(2.*(A10-A11))+
     * 19.*SIN(6.*A10)-9.*SIN(2.*(A10+A10+A11))+7.*COS(2*A10-5.*A13))/
     * 206264.8
	SZH = (24134.*SIN(A11)+5180.*SIN(A10-A11)+4910.*SIN(A10+A11)+
     * 1124.*SIN(A10+A10+A11)+271.*SIN(3.*A10+A11)
     * +132.*SIN(A10+A10-A11)+67.*SIN(4.*A10+A11)+18.*SIN(3.*A10-A11)+
     * 17.*SIN(5.*A10+A11)-10.*SIN(3.*A11)-9.*SIN(A10-3.*A11))/
     * 206264.8
	RH = .39528 -.07834*COS(A10)-.00795*COS(2.*A10)-
     * .00121*COS(3.*A10)-22.E-5*COS(4.*A10)
	CALL HELGEO(DLH,SZH,RH,RA,DEC,DLE,SZE,R)
	 END

	SUBROUTINE WENUS(RA,DEC,DLE,SZE,R)
	COMMON /PLAN/ T
	DOUBLE PRECISION T,D1,D2
	F(D1,D2) = 6.28318530718D0*DMOD(D1 + T*D2,1D0)
	A8  = F(.993126D0,.27377785D-2)
	A13 = F(.140023D0,4.45036173D-3)
	A12 = F(.505498D0,4.45046867D-3)
	A14 = F(.292498D0,4.45040017D-3)
	DLH = A12+(SIN(A13)*(2794.-SNGL(T)/1826.25)-181.*SIN(A14+A14)
     *+12.*SIN(A13+A13) - 10.*COS(2.*(A8-A13)) + 7.*COS(3.*(A8-A13)))
     * /206264.8
	SZH = (12215.*SIN(A14) + 166.*SIN(A13)*COS(A14))/206264.8
	RH  = .72335 - .00493*COS(A13)
	CALL HELGEO(DLH,SZH,RH,RA,DEC,DLE,SZE,R)
	END

	SUBROUTINE MARS(RA,DEC,DLE,SZE,R)                              
	COMMON /PLAN/ T
	DOUBLE PRECISION T,D1,D2
	F(D1,D2) = 6.28318530718D0*DMOD(D1 + T*D2,1D0)
	A8  = F(.993126D0,.27377785D-2)
	A13 = F(.140023D0,4.45036173D-3)
	A15 = F(.987353D0,.145575328D-2)
	A16 = F(.053856D0,.145561327D-2)
	A17 = F(.849694D0,.145569465D-2)
	A19 = F(.056531D0,.23080893D-3)
	DLH = A15+((38451.+(37.+8.*COS(A16))*(1.+SNGL(T)/36525.))*
     * SIN(A16)+2238.*SIN(A16+A16)+181.*SIN(3.*A16)-52.*SIN(2.*A17)-
     * 22.*COS(A16-2.*A19)-19.*SIN(A16-A19)+17.*(COS(A16-A19)+SIN(4.*
     * A16))-16.*COS(2.*(A16-A19))+13.*COS(A8-A16-A16)-10.*(SIN(A16-
     * 2.*A17)+SIN(A16+2.*A17))+7.*(COS(A8-A16)-COS(2.*A8-4.*.75*A16))-
     * 5.*(SIN(A13-3.*A16)+SIN(A8-A16)+SIN(A8-2.*A16))+4.*(-COS((2.*A8-
     * 4.*A16))+COS(A19))+3.*(COS(A13-3.*A16)+SIN(2.*(A16-A19))))
     * /206264.8
	SZH = (6603.*SIN(A17)+622.*SIN(A16-A17)+615.*SIN(A16+A17)+
     * 64.*SIN(2.*A16+A17))/206264.8
	RH = 1.53031-.1417*COS(A16)-.0066*COS(2.*A16)-.00047*COS(3.*A16)
	CALL HELGEO(DLH,SZH,RH,RA,DEC,DLE,SZE,R)
	 END

	SUBROUTINE JOWISZ(RA,DEC,DLE,SZE,R)                             
	COMMON /PLAN/ T
	DOUBLE PRECISION T,D1,D2
	F(D1,D2) = 6.28318530718D0*DMOD(D1 + T*D2,1D0)
	W = T/36525.+1.
	A19 = F(.056531D0,.23080893D-3)
	A25 = F(.400589D0,.3269438D-4)
	A22 = F(.882987D0,.9294371D-4)
	A18 = A19+.2078289
	A0 = 2.*A19-5.*A22
	S = SIN(A19)
	C = COS(A19)
	C0 = COS(A0)
	S0 = SIN(A0)
	S1 = SIN(A19-A22)
	C1 = COS(A19-A22)
	S2 = SIN(A22)
	C2 = COS(A22)
	DLH = A18+(W*(5023.-74*C+68.*S-43.*S0-19.*C0-4.*(C*C-S*S)+6.*S*C
     * -2.*S0*C)+19934.*S+2511.+1093.*C0+1202.*S*C-479*S0-370.*S1*C1+
     * C*(137.*S0-262.*S1*C1)+S*(137.*C0+131.*(C1*C1-S1*S1))+79.*C1-
     * 76.*(C1*C1-S1*S1)+66.*COS(2.*A19-3.*A22)+116.*C0*C-10.*S0*S-37.*C
     * +49.*SIN(2.*A19-3.*A22)+25.*(SIN(A18+A18)+SIN(3.*A19))-
     * 23.*SIN(A0-A19)+17.*(COS(A0+A22)+COS(3.*(A19-A22)))-14.*S1-
     * 13.*SIN(A0+A19+A22)+9.*(C2-COS(A18+A18)-S2-SIN(3.*A19-A22-A22)+
     * SIN(A0+A19+A19)+SIN(2.*(A19-3.*A22)+3.*A25))-8.*(C0*C0-S0*S0)+
     * 7.*(COS(A0+A19+A22)-COS(A19-3.*A22)-2.*S0*C0-SIN(A19-3.*A22))+
     * 6.*(COS(A0+A19+A19)-SIN(3.*(A19-A22)))+5.*(C2*C2-S2*S2)-
     * 4.*(4.*S1*C1*(C1*C1-S1*S1)-C*SIN(A18+A18)+COS(3.*A22)-COS(2.*A19-
     * A22)+COS(3.*A19-2.*A22))+3.*(COS(5.*A22)+COS(5.*(A19-2.*A22))+
     * S2))/206264.8
	SZH = (-4692.*C+259.*S+454.*S*S+W*(30.*S+21.*C)+29.*S*C0+
     * C*(3.*(S0-4.+8.*S*(S+S+1.))+C0+C0)-12.*S*S0)/206264.8
	RH = 5.20883-.25122*C-.00604*(C*C-S*S)+.0026*COS(2.*(A19-A22))+
     * 1.E-5*(-170.*COS(A0+A19)-106.*SIN(2.*(A19-A22))-W*(91.*S+84.*C)+
     * 69.*SIN(A0+A22+A22)-67.*SIN(A0-A19)+63.*S1-51.*COS(A0+A22+A22)-
     * 46.*S+S0*(66.*C-29.*S)+C0*(66.*S-29.*C)+27.*COS(A19-2.*A22)-
     * 22.*C*(1.-4.*S*S)-21.*S0)
	CALL HELGEO(DLH,SZH,RH,RA,DEC,DLE,SZE,R)
	 END

	SUBROUTINE SATURN(RA,DEC,DLE,SZE,R)
	COMMON /PLAN/ T
	DOUBLE PRECISION T,D1,D2
	F(D1,D2) = 6.28318530718D0*DMOD(D1 + T*D2,1D0)
	W = T/36525.+1.
	A19 = F(.056531D0,.23080893D-3)
	A25 = F(.400589D0,.3269438D-4)
	A22 = F(.882987D0,.9294371D-4)
	A21 = A22+1.5727316
	A0 = 2.*A19-5.*A22
	C0 = COS(A0)
	S0 = SIN(A0)
	S = SIN(A22)
	C = COS(A22)
	DLH = A21+(W*(5014.-229.*C-142.*S+101.*S0+60.*C0+22.*SIN(A0+A22)
     * -22.*S*C-17.*(C*C-S*S)-6.*(2.*S0*(S0+C0)-1.-SIN(A0-A22))+
     * 4.*COS(A19-A22-A22)+3.*(COS(2.*A19-4.*A22)-SIN(A21+A21)-COS(A21+
     * A21)+COS(A19+A19-6.*A22))+2.*SIN(A19-A22-A22))+
     *  23045.*S-2689.*C0+2507.-826.*COS(A0+A22)+1604.*S*C+1177.*S0+
     * 425.*SIN(A19-A22-A22)-153.*COS(A0-A22)-114.*C-70.*COS(A21+A21)+
     * 67.*SIN(A21+A21)+66.*SIN(A0-A22)+41.*SIN(A19-3.*A22)+39.*SIN(3.*
     * A22)+31.*(SIN(A19-A22)+ SIN(2.*(A19-A22)))-29.*COS(A19+A19-3.*
     * A22)-28.*SIN(A0-A22+3.*A25)+28.*COS(A19-3.*A22)-
     * 22.*SIN(A22-3.*A25)+20.*(SIN(2.*A19-3.*A22)+COS(2.*A0))
     * +19.*(COS(2.*A22-3.*A25)+2.*S0*C0)-16.*COS(A22-3.*A25)-
     * 12.*(SIN(A0+A22)-COS(A19)+SIN(2.*(A22-A25)))-11.*COS(A0-A22-A22)+
     * 10.*(SIN(A22+A22-3.*A25)+COS(2.*(A19-A22)))+9.*SIN(4.*A19-9.*A22)
     * +8.*(COS(A21*2.-A22)-COS(A21*2.+A22)-SIN(A22-2.*A25)-SIN(2.*A21-
     * A22)+COS(A22+A25))+7.*(SIN(2.*A21+A22)-COS(A19-2.*A22)-COS(A22*2.
     * ))+5.*(SIN(A0-A22-A22)+SIN(3.*A19-4.*A22)-SIN(3.*A19-7.*A22)-
     *COS(3.*(A19-A22))-COS(2.*(A22-A25)))+4.*(SIN(3.*(A19-A22))+SIN(A0+
     * A19))+3.*(COS(2.*A19-6.*A22+3.*A25)+COS(3.*A19-7.*A22)+COS(4.*A19
     * -9.*A22)+SIN(3.*(A19-A22-A22))+SIN(2.*A19-A22)+SIN(A19-4.*A22))
     * +2.*(COS(3.*(A22-A25))+SIN(4.*A22)-COS(3.*A19-4.*A22)-COS(2.*
     * A19-A22)-SIN(2.*A19-7.*A22+3.*A25)+COS(A19-4.*A22)+COS(2.*A0-A22)
     * -SIN(A22-A25)))/206264.8
	SZH = (8297.*S-3346.*C+924.*S*C-189.*(C*C-S*S)+185.+
     * W*(79.*C+18.*S-10.-8.*S*C+3.*SIN(A0+A22))-71.*COS(A0+A22)+
     * 46.*SIN(A0-A22)-45.*COS(A0-A22)+29.*SIN(3.*A22)-
     * 20.*COS(A0+A22+A22)-14.*COS(A0)-11.*COS(3.*A22)+9.*SIN(A19-3.*
     * A22)+8.*SIN(A19-A22)-6.*SIN(A0+A22+A22)+5.*(SIN(A0-A22-A22)-
     * COS(A0-A22-A22))+4.*S0+3.*SIN(A19-A22-A22)*(1.+S+S)+
     * 2.*(SIN(4.*A22)-COS(2.*(A19-A22))))/206264.8
	RH = 9.55774+W*(-.00524*S+.00328*C-.00028)+
     * 1.E-5*(-53252.*C-1878.*SIN(A0+A22)-1482.*(C*C-S*S)+817.*
     * SIN(A19-A22)-539.*COS(A19-A22-A22)+349.*S0+347.*SIN(A0-A22)-
     * 225.*S+149.*COS(A0-A22)-126.*COS(2.*(A19-A22))+104.*COS(A19-A22)+
     * 101.*C0+98.*COS(A19-3.*A22)-73.*COS(A0+A22+A22)-62.*COS(3.*A22)+
     * 42.*SIN(2.*A22-3.*A25)+41.*SIN(2.*(A19-A22))-40.*(SIN(A19-3.*
     * A22)-COS(A0+A22))-23.*SIN(A19)+20.*SIN(A0-A22-A22))
	CALL HELGEO(DLH,SZH,RH,RA,DEC,DLE,SZE,R)
	END

	SUBROUTINE URAN(RA,DEC,DLE,SZE,R)                              
	COMMON /PLAN/ T
	DOUBLE PRECISION T,D1,D2
	F(D1,D2) = 6.28318530718D0*DMOD(D1 + T*D2,1D0)
	A19 = F(.056531D0,.23080893D-3)
	A22 = F(.882987D0,.9294371D-4)
	A25 = F(.400589D0,.3269438D-4)
	A26 = F(.664614D0,.3265562D-4)
	A28 = F(.725368D0,.1672092D-4)
	W = T/36525.+1.
	S = SIN(A25)
	C = COS(A25)
	C2 = COS(A22-A25-A25)
	S2 = SIN(A22-A25-A25)
	S3 = SIN(A22-3.*A25)
	C3 = COS(A22-3.*A25)
	S1 = SIN(A22-A25)
	DLH = A25+2.950458+(19397.*S+570.*2.*S*C+143.*S2+102.*S3+76.*C3-
     * 49.*SIN(A19-A25)+29.*(SIN(2.*A19-6.*A22+3.*A25)+COS(2.*(A25-A28))
     * )-28.*COS(A25-A28)+23.*SIN(3.*A25)-21.*COS(A19-A25)+
     * 20.*(SIN(A25-A28)+C2)-19.*COS(A22-A25)+17.*SIN(2.*A25-3.*A28)+
     * 14.*SIN(3.*(A25-A28))+13.*S1-12.*C+10.*SIN(2.*(A25-A28))-
     * 9.*(SIN(A26+A26)-COS(A25+A25-3.*A28))+8.4853*SIN(2.*(A19-3.*A22+
     * A25+.3926991))+5.*SIN(A22-4.*A25)-4.*(SIN(3.*A25-4.*A28)-
     * COS(3.*(A25-A28)))-3.*COS(A28)-2.*SIN(A28)+
     * W*(110.*S-536*C-30.*(C*C-S*S)+8.*C2+7.*(C3-S3+2.*S*C)+
     * W*(32.-12.*C-9.*S)))/206264.8
	SZH = ((2775.-C)*SIN(A26)+261.*S*COS(A26))/206264.8
	RH = 19.21216-.90154*C-.02121*(C*C-S*S)-.00585*C2 -
     * .00451*COS(A19-A25)+.00336*S1+.00198*SIN(A19-A25)+
     * .00118*C3+.00107*S2-.00081*COS(3.*(A25-A28))-
     * W*(S*.02488+C*.00508+.00206*S*C)
	CALL HELGEO(DLH,SZH,RH,RA,DEC,DLE,SZE,R)
	 END

	SUBROUTINE NEPTUN(RA,DEC,DLE,SZE,R)
	COMMON /PLAN/ T
	DOUBLE PRECISION T,D1,D2
	F(D1,D2) = 6.28318530718D0*DMOD(D1 + T*D2,1D0)
	A19 = F(.056531D0,.23080893D-3)
	A22 = F(.882987D0,.9294371D-4)
	A25 = F(.400589D0,.3269438D-4)
	A28 = F(.725368D0,.1672092D-4)
	A29 = F(.480856D0,.1663715D-4)
	W = T/36525.+1.
	S = SIN(A28)
	C = COS(A28)
	S9 = SIN(A19-A28)
	C9 = COS(A19-A28)
	DLH = (3523.*S-W*(43.*C+4.*S*(1.-W))-50.*SIN(A29+A29)+
     * 29.*S9+38.*S*C-18.*C9+18.38478*SIN(A22-A28+.7853982)-
     * 9.*(SIN(A25+A25-3.*A28)-COS(2.*(A25-A28)))-5.*COS(A25+A25-3.*
     * A28)+4.*COS(A25-A28-A28))/206264.8+A28+.7636834
	SZH = ((6404.-33.*W)*SIN(A29)+110.*S*COS(A29))/206264.8
	RH = 30.07175-W*.00314*S-.25701*C-.00787*COS(A25-A28-A28+4.3735793
     * )+.00409*C9+.0025*S9-.00194*SIN(A22-A28)+.00185*COS(A22-A28)
	CALL HELGEO(DLH,SZH,RH,RA,DEC,DLE,SZE,R)
	 END

	SUBROUTINE PLUTON(RA,DEC,DLE,SZE,R)
	COMMON /PLAN/ T
	DOUBLE PRECISION T,D1,D2
	F(D1,D2) = 6.28318530718D0*DMOD(D1 + T*D2,1D0)
	W = T/36525.+1.
	A31 = F(.663854D0,1.115482D-5)
	A32 = F(.04102D0,1.104864D-5)
	A33 = A32+1.987591
	S = SIN(A32)
	C = COS(A32)
	S1 = SIN(A33)
	C1 = COS(A33)
	DLH = A31+(S*(106892.+S*(S1*C1*(3900.+2280.*C)-S*(C*6712+12516.))
     * +C*33312.+W*(W*227.+200.)+S1*S1*(7574.+C*(2156.+2280.*C)))-
     * S1*C1*(9136.-  90.*C))*4.8481368E-6
	SZH = (57726.*S1+29359.*S*C1-1155.*C*S1+3870.*SIN(A32+A32+A33)+
     * 1138.*SIN(3.*A32+A33)+472.*SIN(A32+A32-A33)+353.*SIN(4.*A32+A33)-
     * 255.*S*C1*(1.-4.*S1*S1)+(33.*C-119.)*S1*(4.*C1*C1-1.))/206264.8
	RH = 41.86345-C*(8.90288+C*(1.93438+C*(.90596+.39968*C)))
	CALL HELGEO(DLH,SZH,RH,RA,DEC,DLE,SZE,R)
	 END

	subroutine H

C   EFEMERYDY UKLADU SLONECZNEGO PLUS HOROSKOP

	DIMENSION RA(10),D(10),DL(10),B(10),R(10)
	character*1 KOMENT*80,NAZWA(9),YES
	write(*,1)
1	FORMAT(1H0/10X,'PROGRAM "HORR"'//
     *' OBLICZA EFEMERYDY PLANET, ICH ASPEKTY DLUGOSCIOWE I ',
     *'DEKLINACYJNE ORAZ ROZKLAD'/' NA ZNAKI ZODIAKU I DOMY ',
     *'ASTROLOGICZNE. PONADTO PODAJE TABELE ZYWIOLOW I JAKOS-'
     */' CI, ILOSC PLANET POZYTYWNYCH I NEGATYWNYCH (DOLNY PRAWY'
     *,' ROG TEJ TABELI) ORAZ'/' ILOSCI PLANET WIDOCZNYCH I NIE'
     *,'WIDOCZNYCH (NAD I POD HORYZONTEM)'/
     *' I PLANET WSCHODNICH I ZACHODNICH (NA KONCU ROZKLADOW PLANET'
     *,' NA DOMY).'/
     *'   DWA PODZIALY NA DOMY PRZYJETO WEDLUG PLACIDUSA (PIERWSZY)'
     *,' I C. EVANSA (TZW.'/' NATURALNEJ GRADACJI; DRUGI PODZIAL).')
	write(*,2)
2	FORMAT(1H /'    WYMAGANE DANE WEJSCIOWE:'/' (1) DATA (ROK,'
     *,' MIESIAC, DZIEN, GODZINA UT)'/' (2) WSPOLRZEDNE: DLUGOSC'
     *,' (>0 NA WSCHOD) I SZEROKOSC GEOGR. (STOPNIE)'/
     *' (3) DO 80 ZNAKOW TEKSTU, z ktorego pierwsze 9 znakow',
     *' bedzie BLOKOWYM tytulem'/)
	data ifirst/0/
	if(ifirst.eq.0) open(6,file='h.dat',status='UNKNOWN')
	ifirst=ifirst+1
10	write(*,3)
3	FORMAT(1H ,'   Podaj teraz swoje dane (zaraz po tym przyk'
     *,'ladzie)'/' 1926,11,23,2359,+18.5598,+53.098'/' ANNAMARIA'
     *,' ABACKA, PIWNICE K. TORUNIA, CESARSKIE CIECIE'/)
	read(*,*) IR,M,ID,IUT,DLG,SZG
	READ(*,'(A80)') KOMENT
	READ(KOMENT,'(9A1)') NAZWA
	write(*,5) IR,M,ID,IUT,DLG,SZG,KOMENT
5	FORMAT(1H ,'   WCZYTANO NASTEPUJACE DANE:'/4I5,2F10.5/
     *1X,a/' O.K.? [ENTER=tak]: ')
	read(*,'(a)') YES
	IF(YES.nE.' '.and.YES.NE.'Y'.AND.YES.NE.'y') GO TO 10
	call fpage(NAZWA,9)
	WRITE(6,14)      ID,M,IR,IUT,DLG,SZG,KOMENT
	UT=IUT/100+MOD(IUT,100)/60.
14	FORMAT(1H //10X,'DANE: ',2I3,I5,I6,2F10.4,10X,A/)
	DLGr=DLG*.1745329E-1
	SZGr=SZG*.1745329E-1
	CALL EFEMER(IR,M,ID,UT,RA,D,DL,B,R)
	CALL HOR(DLGr,SZGr,DL,D)
	WRITE(6,15)
	write(*,16)
16	format(1h ,7x,' Wyniki upchnalem w zbiorze ===> H.DAT'/ 
     *' chcesz kontynuowac obliczyc nast. horoskop? [ENTER=nie] ')
	read(*,'(A)') YES
	IF(YES.nE.' '.and.YES.ne.'Y'.and.YES.nE.'y') GO TO 10
	close(6)
15	FORMAT(1H /)
	END

	SUBROUTINE EFEMER(IY,M,ID,UT,RA,D,DL,B,R)
	DOUBLE PRECISION DJ
	DIMENSION RA(10),D(10),DL(10),B(10),R(10)
	DMS(Y)=AINT(Y)+AINT(60.*AMOD(Y,1.))*.01+AMOD(60.*Y,1.)*.006
c        JD(L,M,N) = (L+(M-8)/6)*1461/4+(153*MOD(M+9,12)+2)/5+N+1721119
c     *                                        -(L+100+(M-8)/6)/100*3/4
	B(1)=0.
	j=jd(iy,m,id,0)
	DJ=j-.5D0+UT/24.
	CALL Sun(dj,RA(1),D(1),DL(1),R(1),RC,st)
	CALL KSIEZY(dj-.01,RA1,D1,DL1,B1,R1)
	CALL KSIEZY(dj+.01,RA1,D1,DL2,B2,R1)
	CALL KSIEZY(dj,RA(2),D(2),DL(2),B(2),R(2))
	OM=SIN(B(2))/COS(B(2))
	OM=ATAN2(OM,SIGN(SQRT(abs(8.122931E-3-OM*OM)),B2-B1))
	OM=AMOD((DL(2)-OM)*57.29578+360.,360.)
	WRITE(6,9) 
	II=MOD(J+1,7)
	WRITE(6,10) ID,M,IY,J,UT,OM,II
	DO 2 I=3,10
2	CALL PLANET(I-2,RA(I),D(I),dl(I),b(i),r(i))
	 DO 1 I=1,10
	RA1=DMS(RA(I)*3.8197186)
	D1=DMS(D(I)*57.295779)
	DL1=DMS(DL(I)*57.295779)
	B1 = DMS(B(I)*57.295779)
1	 WRITE(6,11) MOD(I,10),RA1,D1,DL1,B1,R(I)
	RETURN
9	FORMAT(1H //,'     EFEMERYDY CIAL UKLADU SLONECZNEGO (JD JEST DATA
     * JULIANSKA W POLUDNIE UNIWERSALNE, TZN. O GODZ. 12 UT)',/,4X,
     *'W KOLEJNYCH KOLUMNACH PODANO WSPOLRZEDNE ROWNIKOWE I EKLIPTYCZNE
     * (GEOCENTRYCZNE) ORAZ ODLEGLOSC OD ZIEMI DO: SLONCA (NR. 1), ',
     */,4X,'KSIEZYCA (2), MERKUREGO (3) I T.D. AZ DO PLUTONA (0)',/)
10	FORMAT( /'   DNIA:',2I3,I5,'   JD =',I8,'  GODZ.',F6.3,' UT',
     *T70,' DLUGOSC WEZLA WSTEPUJACEGO KSIEZYCA =',F7.2,' STOPNI'/
     *'   DZIEN TYGODNIA =',I2,' (TUTAJ 0 = NIEDZIELA, 1 = PONIEDZ..)'//
     *'  NR    RA       DEKL.     DLUG.     SZER.     ODLEG.',/,
     *'     GODZ.MMSS STOP.MMSS STOP.MMSS STOP.MMSS   AU/RZ')
11	FORMAT( I4,5F10.4,3F15.5,I4,A8)
	END

	SUBROUTINE HOR(DLG,SZERG,DL,D)                             
	COMMON /suncom/ xi0,eta0,SE,CE,dla /plan/ dj
	character*8 pl(10),ajak(3),sp*1,zer*1,plu*1
	DOUBLE PRECISION dj,T
	diMENSION C(12),L(12),DL(10),D(10),plu(12)
	DATA PL/' SLONCE ',' KSIEZYC',' MERKURY',' WENUS  ',' MARS   ',
     *' JOWISZ ',' SATURN ',' URAN   ',' NEPTUN ',' PLUTON '/
     *,AJAK/'KARDYNAL','STALE   ','ZMIENNE '/,sp,zer/' ','0'/
	PI=3.1415927
	DG=180./PI
	T=dj/36525D0
	ST=AMOD((SNGL(DMOD(673.1054841D0/864D0+dj+T/864D2*
     *(8640184.812866D0+T*(.093104D0-6.2D0*T)),1D0))+2.)*(PI+PI)+DLG,
     *PI+PI)
	WRITE(6,9) PL
	TD=AMOD(PI+PI+ATAN2(COS(ST),-SE*TAN(SZERG)-CE*SIN(ST)),PI+PI)
	TM=AMOD(PI+PI+ATAN2(SIN(ST),COS(ST)*CE),PI+PI)
	DO 3 I=1,10
	II=I-1
	DO 4 J=1,II
4	 C(J)=(D(J)-D(I))*DG
	DO 2 J=I,12
	IF(J.LT.11) A=DL(J)
	IF(J.EQ.11) A=TD
	IF(J.EQ.12) A=TM
	A=ABS(A-DL(I))*DG
	IF(A.GT.180.) A=360.-A
2	C(J)=A
3	WRITE(6,1) MOD(I,10),PL(I),C
	WRITE(6,21) ((I-1)*30,I=1,12)
	DO 6 I=1,12
	IF(I.LT.11) DL(I)=DL(I)*DG
	C(I)=0.
	PLU(I)=SP
6	L(I)=0
	DO 16 I=1,10
	LL=IFIX(DL(I)/30.)+1
	C(LL)=C(LL)+1.
16	IF(I.LT.10) L(LL)=L(LL)*10+I
c	L(LL)=9-L(LL)
	plu(ll)=zer
	WRITE(6,22) (L(i),plu(i),i=1,12),(MOD(I,10),PL(I),I=1,10),
     *DL,TD*DG,TM*DG
	WRITE(6,29)
	DO 11 I=1,9,4
	A=0.
	DO 110 K=1,10,3
110	A=A+C(K+I/4)
11	WRITE(6,28) AJAK(I/4+1),(L(MOD(I+J,12)+1),PLU(MOD(I+J,12)+1),
     * J=35,62,9),INT(A)
	LL=0
	DO 12 J=1,4
	L(J)=0
	DO 12 I=J,12,4
	IF(MOD(J,2).EQ.1) LL=LL+C(I)
12	L(J)=L(J)+C(I)
	K=LL-10
	WRITE(6,30) (L(I),I=1,4),LL,K
	TD=AMOD(PI+PI+ST-ATAN2(-CE*SIN(TD),-COS(TD)),PI+PI)
	DO 5 I=7,12
	C(I)=ABS(FLOAT(10-I))*(TD-FLOAT(I/10)*PI)/3.
	C(I-6)=ATAN2(SIN(ST-C(I))/CE,COS(ST-C(I)))*DG+180.
5	C(I)=AMOD(C(I-6)+180.,360.)
	WRITE(6,26) ST*3.81972,DLG*DG,SZERG*180./PI
	DO 17 K=1,2
	DO 14 I=1,12
	PLU(I)=SP
14	L(I)=0
	LL=0
	DO 7 J=1,12
	A=C(MOD(J,12)+1)
	DO 7 I=1,10
	IF(C(J).GT.A.AND.(DL(I).GE.A.AND.DL(I).LT.C(J))) GO TO 7
	IF(C(J).LT.A.AND.(DL(I).GE.A.OR.DL(I).LT.C(J))) GO TO 7
	IF(I.LT.10) L(J)=L(J)*10+I
	IF(I.EQ.10) PLU(J)=ZER
	IF(K.EQ.1.AND.J.GT.6) LL=LL+1
	IF(K.EQ.2.AND.(J.LT.4.OR.J.GT.9)) LL=LL+1
7	CONTINUE
	J=10-LL
	if(k.eq.1) then
	WRITE(6,25) 'Placidusa: ',C,LL,(L(I),PLU(I),I=1,12),J
c	IF(K.EQ.2) GO TO 17
	A=AMOD(360.+C(1)-C(10),360.)
	DO 180 I=2,6
	C(I)=AMOD(C(1)+(11+I*(I-8))*A/5.-I*(18*I-174)-228,360.)
180	C(I+6)=AMOD(C(I)+180.,360.)
	else
	WRITE(6,25) 'C. Evansa: ',C,LL,(L(I),PLU(I),I=1,12),J
	endif
17	CONTINUE
	 RETURN
1	FORMAT(I4,A8,12F10.2)
9	FORMAT(1H /'   DEKLINACYJNE I DLUGOSCIOWE ASPEKTY SLONCA, KSIEZYCA
     *, PLANET, ASCENDENTU I MEDIUM COELI (W STOPNIACH)',//,
     *'  NR OBIEKT: ',10(2X,A8),' ASCENDENT MEDIUM C.')
21	FORMAT(1H /'   ROZKLAD PLANET NA ZNAKI ZODIAKU',//,4X,5HZNAK:,6X,'
     *BARAN       BYK BLIZNIETA       RAK       LEW     PANNA',
     *6X,'WAGA  SKORPION  STRZELEC KOZIOROZEC   WODNIK      RYBY',/,
     *' DL.EKL.: ',12I10)
22	FORMAT(1H ,'PLANETY: ',12(I9.0,a1)//'   GDZIE:',10(I2,A8),
     *'  (ASC.     M.C. )',/,' DL.EKL.:',F7.2,9F10.2,4X,2F8.2)
26	FORMAT(///,4X,'ROZKLAD PLANET NA DOMY (CZAS GWIAZDOWY =',F7.3,
     *' GODZ., WSPOLRZEDNE GEOGR.:',F7.3,2HE,,F8.3,'N (STOPNI))',/)
25	FORMAT(1H /' DOM ',a11,1HI,8X,2HII,7X,3HIII,8X,2HIV,9X,1HV,8X,
     *2HVI,7X,3HVII,6X,4HVIII,8X,2HIX,9X,1HX,8X,2HXI,7X,3HXII,/,
     *9H DL.EKL.:,12F10.2,I3,/,9H PLANETY:,12(I9.0,A1),I3)
28	FORMAT(5H   I ,A8,4(2H *,I9.0,A1),2H =,I3,4H   I)
29	FORMAT( 1H /'   I JAKOSC|TRYGON: OGIEN +    ZIEMIA + POWIETRZE +'
     *,'      WODA = SUMA I')
30	FORMAT(13H   I    RAZEM ,4(2H &,I10),2H =,I2,I3,2H I)
      END



	SUBROUTINE FPAGE(text,N)

C     PRINTS OUT A LINE OF SYMBOLS BUILT UP OF THE LINEPRINTER
C     CHARACTERS. A SINGLE LINE CAN CONTAIN MAX. 10 SYMBOLS. THE NUMBERS
C     FROM 0 TO 9, THE FULL ALPHABET, AND THE SYMBOLS '.', ':', '-', ' '
C     ARE AVAILABLE. AN ARBITRARY TEXT CAN BE PRINTED OUT BY REPEATABLE
C     CALLS TO FPAGE.

C     text - ARRAY OF SYMBOLS TO BE PRINTED OUT IN ONE LINE. THERE CAN
C            BE MAX. 10 CHARACTERS IN text (I),
C     N -    NO. OF SYMBOLS IN text.

	character*1 BYTE,ISYMB(40),text(10),blok*5280
	DIMENSION ISEL(10)
	DATA ISYMB/1HA,1HB,1HC,1HD,1HE,1HF,1HG,1HH,1HI,1HJ,1HK,1HL,1HM,
     *1HN,1HO,1HP,1HQ,1HR,1HS,1HT,1HU,1HV,1HW,1HX,1HY,1HZ,
     *1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9,1H0,1H:,1H.,1H-,1H /
       Blok(1189:2376)=                        '########################
     *         ###         ###         ###         ###         ###      
     *   ######      ############### ########## ###      ######     ####
     *###    #### ###   ####  #########   ########    #########   ###   
     *####  ###    #### ###     #######      ######         ###         
     *###         ###         ###         ###         ###         ###   
     *      ###         ###########################      #######    ####
     *#####  #################### #### ######  ##  ######      ######   
     *   ######      ######      ######      ######      #######     ###
     *#####    #########   ##########  ###### #### ######  ##########   
     *#########    ########     #######      ### ########## ############
     *###      ######      ######      ######      ######      ######   
     *   ######      ############### ########## ########### ############
     *###      ######      ######      ########################## ###   
     *      ###         ###         ###          ########## ############
     *###      ######      ######      ######      ######      ######   
     *## ######    ################  ######## ############# ############
     *###      ######      ######      ########################## ###   
     *###   ###    ###  ###     ### ###      ###'
       Blok(2377:3564)=                        ' ########## ############
     *###         ###         ###########  ###########         ###      
     *   ###         ############### ########## ########################
     *    ####        ####        ####        ####        ####        ##
     *##        ####        ####        ####    ###      ######      ###
     *###      ######      ######      ######      ######      ######   
     *   ######      ############### ########## ###      ######      ###
     *###      ######      ######      ######      ### ###    ###   ### 
     * ###     ######       ####         ##     ###      ######      ###
     *###      ######      ######      ######  ##  ###### #### #########
     *###########  #########    #######      ######      ######      ###
     * ###    ###   ###  ###     ######       ####       ######     ### 
     * ###   ###    ### ###      ######      ######      ######      ###
     *###      ######      ### ###    ###   ###  ###     ######       ##
     *##        ####        ####        ####    ########################
     *       ####       ####       ####       ####       ####       ####
     *       ####       ########################      88         888    
     *    8888       88888      888888        8888        8888        88
     *88        8888      88888888    88888888  '
           Blok(3565:4752)                =    ' 8888888888 888888888888
     *888     8888888    8888       8888       8888       8888       888
     *8       8888       88888888888888888888888 8888888888 888888888888
     *888      888888      888      888888     888888       888888888   
     *   888888      888888888888888 8888888888     888888     8888888  
     *  888  888   888   888  888    888  888888888888888888888888      
     * 888         888         888         888  888888888888888888888888
     *888         888         88888888888 888888888888         888      
     *   888888      888888888888888 8888888888  8888888888 888888888888
     *888         888         88888888888 888888888888888      888888   
     *   888888      888888888888888 8888888888 888888888888888888888888
     *        888        888        888        888        888        888
     *        888        888        888          8888888888 888888888888
     *888      888888      888888888888888 8888888888 888888888888888   
     *   888888      888888888888888 8888888888  8888888888 888888888888
     *888      888888      888888      888888888888888 88888888888      
     *   888         888888888888888 8888888888   88888888   8888888888 
     *888     8888888    88888888   88 888888  88  888888 88   88888888 
     *   8888888     888 8888888888   88888888  '
          Blok(4753:5280) =                    '                        x
     *    ####        ####        ####                                ##
     *##        ####        ####                                        
     *                                                                  
     *          ####        ####        ####                            
     *                                     ##########  ##########       
     *                                                                  
     *                                                                  
     *                                          '
        Blok(1:1188)=                          ' ########## ############
     *###      ######      ######      ##############################   
     *   ######      ######      ######      ############## ############
     *###      ######      ############## ##########  ########### ###   
     *   ######      ##########################  ########## ############
     *###      ######      ######         ###         ###         ###   
     *   ######      ############### ########## ##########  ########### 
     *###     #######      ######      ######      ######      ######   
     *   ######     ############### ##########  ########################
     *###         ###         ###         #########   #########   ###   
     *      ###         ################################################
     *###         ###         ###         #########   #########   ###   
     *      ###         ###         ###          ########## ############
     *###      ######      ######         ###         ###    ########   
     * ########      ############### ########## ###      ######      ###
     *###      ######      ######      ##############################   
     *   ######      ######      ######      ###########################
     *    ####        ####        ####        ####        ####        ##
     *##        ####    ########################'

	IF(N.LE.0 .OR. N.GT.10) GO TO 10
c	write(6,'(a5280)') blok


C    IDENTIFY THE CHARACTERS (LOWERCASE LETTERS ARE UPPERCASED)

2	DO 8 J=1,N
	BYTE=text(J)
	IF(ichar(BYTE).GT.96.AND.ichar(BYTE).LT.123)
     *                     BYTE=char(ichar(BYTE)-32)
	DO 7 K=1,39
	IF(BYTE.Eq.ISYMB(K)) GO TO 8
7	CONTINUE
8	ISEL(J)=K

C    OUTPUT THE BLOCK LETTERS

	do 20 lw=1,11
	l=1+(lw-1)*12
20	write(6,'(1x,10(a12,3x))') 
     *(blok((ISEL(J)-1)*132+l:(ISEL(J)-1)*132+L+11),J=1,N)

	RETURN

10	WRITE(*,'(A,I3,A)') '*** Wrong value of parameter N (',N,
     *       ') passed to subroutine FPAGE.'
	STOP 16
	END




	Subroutine Rocznik(iy)
c do obliczania dat swiat w danym roku
	character MIES(12)*12,dtyg(7)*4,swieto(50)*38 !ofile*10,
	dimension js(50)
	real*8 jswieto(50),j0
	DATA dtyg/4hNie ,4hPon ,4hWto ,4h/Sro,4hCzw ,4hPi/a,4hSob /
     * mies/'stycze/n',4hluty,8hmarzec  ,'kwiecie/n',3hmaj,
     * 8hczerwiec,8hlipiec  ,'sierpie/n',9hwrzesie/n,'pa/xdziernik',
     * 8hlistopad,'grudzie/n'/
c	iy=2010
c	write(ofile,'(1hR,i4,a)') iy,'.txt'
c	open(1,ofile)

	write(6,'(/a,i4,a)') '\break\ziel\centering'
     &//'\LARGE\bf /Swi/eta pa/nstwowe i ko/scielne w Polsce\\ oraz'
     &//' inne wa/zniejsze dni w roku ',iy,'\\[5mm]'
	write(6,'(a)') '\bigskip\setbox0\vbox{\lineskip=0pt\large\halign'
     &//'{#\ \hfill&#\hfill&\hfil#&#\hfill&\quad\strut\vrule\,\vrule\ '
     &//'#&#\hfill&#\hfil&\hfill#&#\cr'

	NW = NWIELK(iy)
	N=NW/10
	M=mod(NW,10)
c	print*,NW,N,M
	J=JD(iy,M,N,0)
	 swieto(1)='Julia/nski Nowy Rok '
	jswieto(1)=JD(iy,1,1,1)
c	swieto(1)='Siedemdziesi/atnica'
c	jswieto(1)=J-63
	swieto(3)='/Sroda Popielcowa'
	jswieto(3)=J-46
	swieto(2)='Niedziela Palmowa'
	jswieto(2)=J-7
	swieto(9)='Niedziela Wielkanocna'
	jswieto(9)=j
	swieto(4)='Wniebowstapienia Pa/nskiego'
	jswieto(4)=J+39
	swieto(5)='Zes/lania Ducha /Swi/etego'
	jswieto(5)=J+49
	swieto(7)='Tr/ojcy Naj/swi/etszej'
	jswieto(7)=j+56
	swieto(6)='Bo/ze Cia/lo'
	jswieto(6)=J+60
	swieto(8)='Naj/swi/etszego Serca Pana Jezusa'
	jswieto(8)=J+68
	swieto(10)='T/lusty czwartek'
	jswieto(10)=J-52
	swieto(11)='Nowy Rok'
	jswieto(11)=JD(iy,1,1,0)
	swieto(12)='Objawienia Pa/nskiego (Trzech Kr/oli)'
	jswieto(12)=JD(iy,1,6,0)
	swieto(13)='Ofiarowania Pa/nskiego'
	jswieto(13)=JD(iy,2,2,0)
	swieto(14)='Naj/swi/etsza Maryja Kr/olowa Polski'
	jswieto(14)=JD(iy,5,3,0)
	swieto(15)='Wszystkich /Swi/etych'
	jswieto(15)=JD(iy,11,1,0)
	swieto(16)='Zaduszki'
	jswieto(16)=JD(iy,11,2,0)
	swieto(17)='Wniebowzi/ecia N.M.P.'
	jswieto(17)=JD(iy,8,15,0)
	swieto(18)='Bo/ze Narodzenie'
	jswieto(18)=JD(iy,12,25,0)
	swieto(19)='/Swi/eto Pracy'
	jswieto(19)=JD(iy,5,1,0)
	swieto(20)='/Swi/eto Konstytucji 3 Maja'
	jswieto(20)=JD(iy,5,3,0)
	swieto(21)='/Swi/eto Niepodleg/lo/sci'
	jswieto(21)=JD(iy,11,11,0)
	swieto(22)='Naj/swi/etszej Maryi Panny Rodzicielki'
	jswieto(22)=JD(iy,1,1,0)!Uroczystoœæ Œwiêtej Bo¿ej Rodzicielki Maryi
	swieto(23)='Dzie/n Flagi Rzeczpospolitej Polskiej'
	jswieto(23)=JD(iy,5,2,0)
	swieto(24)='/Swi/eto Wojska Polskiego'
	jswieto(24)=JD(iy,8,15,0)
	swieto(25)='Dzie/n Solidarno/sci i Wolno/sci'
	jswieto(25)=JD(iy,8,31,0)
	swieto(33)='Dzie/n Papie/za Jana Paw/la II'
	jswieto(33)=JD(iy,10,16,0)
	it=MOD(int(jswieto(33))+1,7)+1
	 swieto(26)='Dzie/n Papieski' ! niedziela poprzedz. 16 X
	jswieto(26)=jswieto(33)-it+1

	 swieto(27)='Dzie/n Babci'
	jswieto(27)=JD(iy,1,21,0)
	 swieto(28)='Dzie/n Dziadka'
	jswieto(28)=JD(iy,1,22,0)
	 swieto(29)='Matki Boskiej Gromnicznej'
	jswieto(29)=JD(iy,2,2,0)
	 swieto(30)='Dzie/n Matki'
	jswieto(30)=JD(iy,5,26,0)
	 swieto(31)='Dzie/n Ojca'
	jswieto(31)=JD(iy,6,23,0)
	 swieto(32)='Noc /Swi/etoja/nska'
	jswieto(32)=JD(iy,6,23,0)
	 swieto(34)='Pierwsza niedziela Adwentu'
	 it=MOD(int(jswieto(18))+1,7)+1
	 if(it.eq.1) it=it+7
	jswieto(34)=jswieto(18)-it-20
	 swieto(35)='Sylwester'
	jswieto(35)=JD(iy,12,31,0)
	 swieto(36)='Niedziela Mi/losierdzia Bo/zego'
	jswieto(36)=J+7
	 swieto(37)='/Swi/eto UMK'
	jswieto(37)=JD(iy,2,19,0)
	 swieto(38)='Walentynki'
	jswieto(38)=JD(iy,2,14,0)
	 swieto(39)='Dzie/n Kobiet'
	jswieto(39)=JD(iy,3,8,0)
	 swieto(40)='Narodowy Dzie/n /Zycia'
	jswieto(40)=JD(iy,3,24,0)
	 swieto(41)='Dzie/n Edukacji Narodowej'
	jswieto(41)=JD(iy,10,14,0)
	 swieto(42)='Dzie/n Polskiej Niezapominajki'
	jswieto(42)=JD(iy,5,15,0)
	 swieto(43)='Dzie/n Dziecka'
	jswieto(43)=JD(iy,6,1,0)
	 swieto(44)='/Sw.\ Piotra i Paw/la'
	jswieto(44)=JD(iy,6,29,0)
	 swieto(45)='Niepokalane Pocz/ecie N.M.P.'
	jswieto(45)=JD(iy,12,8,0)
	 swieto(46)='Do/zynki - /Swi/eto Plon/ow'
	jswieto(46)=JD(iy,9,23,0)
	 it=mod(7-MOD(JD(iy,6,1,0)+1,7),7)
	 swieto(47)='Dzie/n Dzi/ekczynienia'
	jswieto(47)=JD(iy,6,1,0)+it
	 swieto(48)='Dzie/n Programisty'
	jswieto(48)=JD(iy,9,13,0)

	nsw=48

c Na wypadek, gdy kilka swiat zbiega sie w tym samym dniu
	do 20 i=1,nsw-1
	do 20 j=i+1,nsw
20	if(jswieto(i).eq.jswieto(j)) jswieto(j)=jswieto(j)+.1

c Odszukanie najwczesniejszej daty
	j0=1d7
	do 1 i=1,nsw
	if(jswieto(i).lt.j0) then
		js(1)=i
		j0=jswieto(i)
	endif
1	continue

c Sortowanie wg dat, ustawianie kolejnosci swiat w kalendarzu
	do 2 i=2,nsw 
	j0=1d7
	do 2 j=1,nsw
	if(jswieto(j).lt.j0.and.jswieto(j).gt.jswieto(js(i-1))) then
		js(i)=j
		j0=jswieto(j)
	endif
2	continue

	ns=nsw/2+mod(nsw,2)
	do 10 j=1,ns
	i=js(j)
	it=MOD(int(jswieto(i))+1,7)+1
	call DATA(int(jswieto(i)),ir,im,id,0)
c10	write(6,'(a,1x,a4,i3,1x,a12,f12.1)') swieto(i),dtyg(it),id,
c     & Mies(im) !,jswieto(i)
	write(6,'(a,1h&,a4,1h&,i3,1h&,a9,2h&&)') swieto(i),dtyg(it),id,
     & Mies(im)
		if(j+ns.le.nsw) then
	i=js(j+ns)
	it=MOD(int(jswieto(i))+1,7)+1
	call DATA(int(jswieto(i)),ir,im,id,0)
	write(6,'(a,1h&,a4,1h&,i3,1h&,a12,3h\cr)') swieto(i),dtyg(it),id,
     & Mies(im)
     		else
     	 write(6,'(a)') '&&&\cr'
		endif
10	continue
c	if(mod(nsw,2).eq.1)
	write(6,'(a)') '}}\medskip\centerline{\fbox{\fbox{\box0}}}'

	end
