	program risetnew
	character miasto*12
	dimension ws(12),za(12)
	iTRANS(X)=(2.2918312*X+.4*FLOAT(IFIX(X*3.819719)))*100.
	write(*,'(a,a)') ' Wschody i zachody Slonca dla 28 lat. ',
     *' Podaj nazwe miejsca [a12]: '
	read(*,'(a12)') miasto
	write(*,'(a,a/a)') ' Wspolrzedne geogr. FI, DL [deg], Rok odniesi'
     *,'enia (ujemny rok oznacza BRAK ',
     *'    zapisu wynikow do pliku RISET.DAT):  '
	read(*,*) fid,dld,iy00
	fi=abs(fid)*3.14159265D0/180.
	dl=dld*3.14159265D0/180.
	if(iy00.gt.0) open(6,file='riset.dat',status='unknown')
	iy0=iabs(iy00)
	iy1=iy0+3
	do 200 iy=iy0,iy1
	write(6,'(1h //14x,''Wschody [w] i zachody [z] Slonca'',
     *'' w latach'',i5,6(1h,,i5),''. Miejsce: '',a12/6x,
     *''Dlugosc i szerokosc geogr.:'',2f6.2,'' [stopni]. '',''Momenty ''
     *,''podano w godzinach i minutach czasu srodkowoeuropejskiego''/)')
     * (i-13+iy,i=1,25,4),miasto,dld,fid
	write(6,110) 
110	format(1h ,'    Styczen     Luty     Marzec   Kwiecien    Maj '
     *,'    Czerwiec   Lipiec   Sierpien Wrzesien Pazdziernik'
     *,' Listopad  Grudzien'/2x,'D',2x,12(' w    z   '),'D') 
	do 200 id=1,31
	if(mod(id-1,5).eq.0) write(6,'(1h )')
	do 120 im=1,12
	d0=jd(iy,im,id,0)
	call sunset(d0,dl,fi,1,rc,ws(im),p,za(im),r,d,s,x)
	if(id.lt.29) go to 120
	call data(int(d0),l,m,n,0)
	if(im.eq.m) go to 120
	ws(im)=0.
	za(im)=0.
120	continue
200	write(6,'(1x,i2,12(2i5.0),i3)') id,(itrans(ws(i)+2.181662e-3)
     *,itrans(za(i)+2.181662e-3),i=1,12),id
	if(iy00.gt.0)write(*,'(a)')'        Patrz ============> RISET.dat'
	stop
	end

	SUBROUTINE SUNSET(AJ,DLUGWS,SZEROK,ISTREF,RCZ,
     * WSCH,POLUDN,ZACH,RA,DEKL,RS,DLU)
	common /sun/ xi0,eta0,s,ce,dla
	DOUBLE PRECISION DJD
	DJD=AJ-DBLE(ISTREF/24.)
	CALL earth(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)
	RETURN
	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
	RETURN
	  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
        RETURN
          END

	SUBROUTINE EARTH(DJE,ra,adec,adlu,r,EQT,GST)

C GIVEN  DJE = JULIAN EPHEMERIS DATE.
C OUTPUT REFERRED TO MEAN EQUATOR AND EQUINOX OF DATE (DJE).

C WG: P. STUMPFF, 1980, ASTRON. ASTROPHYS. SUPPL. SER. 41, 1-8.
	IMPLICIT REAL*8 (D)
	DIMENSION  F(3),SN(4)
	DIMENSION DCFEL(3,8),DCEPS(3),CCSEL(3),DCARGS(2,15),
     *	CCAMPS(5,15),CCSEC(3,4),DCARGM(2,3),CCAMPM(4,3)
C dla "sunset"
	common /sun/ xi0,eta0,sep,ce,ala
	DATA DC2PI/6.2831853071796D0/,CC2PI/6.283185/,
     * DC1/1.0D0/,DCT0/2415020.D0/,DCJUL/36525D0/, CCIM /8.978749E-2/

C CONSTANTS DCFEL(I,K) OF FAST CHANGING ELEMENTS
C     DCFEL(I,1)=L srednia dl. Slonca
	DATA DCFEL/1.7400353D+00, 6.2833195099091D+02, 5.2796D-6,
     *		 6.2565836D0,   6.2830194572674D2,  -2.6180D-6,
     *		 4.7199666D0,   8.3997091449254D3,  -1.9780D-5,
     *		 1.9636505D-1,  8.4334662911720D3,  -5.6044D-5,
     *		4.1547339D0,    5.2993466764997D1,   5.8845D-6,
     *		4.6524223D0,    2.1354275911213D1,   5.6797D-6,
     *		4.2620486D0,    7.5025342197656D0,   5.5317D-6,
     *		1.4740694D0,    3.8377331909193D0,   5.6093D-6/
C  sredn. anom. Ks. 5.1680003d0,    1717915922.633"t   31.31"t**2
	DATA DCEPS/4.093198D-1,-2.271110D-4,-2.860401D-8/
	DATA CCSEL/1.675104E-2,-4.179579E-5,-1.260516E-7/

	DATA DCARGS/5.0974222D0,-7.8604195454652D2,
     *3.9584962D0,-5.7533848094674D2,
     *1.6338070D0,-1.1506769618935D3,
     *2.5487111D0,-3.9302097727326D2,
     *4.9255514D0,-5.8849265665348D2,
     *1.3363463D0,-5.5076098609303D2,
     *1.6072053D0,-5.2237501616674D2,
     *1.3629480D0,-1.1790629318198D3,
     *5.5657014D0,-1.0977134971135D3,
     *5.0708205D0,-1.5774000881978D2,
     *3.9318944D0, 5.296346478D1,
     *4.8989497D0, 3.9809289073258D1,
     *1.3097446D0, 7.7540959633708D1,
     *3.5147141D0, 7.9618578146517D1,
     *3.5413158D0,-5.4868336758022D2/

	DATA CCAMPS/
     *-2.279594E-5, 1.407414E-5, 8.273188E-6, 1.340565E-5,-2.490817E-7,
     *-3.494537E-5, 2.860401E-7, 1.289448E-7, 1.627237E-5,-1.823138E-7,
     * 6.593466E-7, 1.322572E-5, 9.258695E-6,-4.674248E-7,-3.646275E-7,
     * 1.140767E-5,-2.049792E-5,-4.747930E-6,-2.638763E-6,-1.245408E-7,
     * 9.516893E-6,-2.748894E-6,-1.319381E-6,-4.549908E-6,-1.864821E-7,
     * 7.310990E-6,-1.924710E-6,-8.772849E-7,-3.334143E-6,-1.745256E-7,
     *-2.603449E-6, 7.359472E-6, 3.168357E-6, 1.119056E-6,-1.655307E-7,
     *-3.228859E-6, 1.308997E-7, 1.013137E-7, 2.403899E-6,-3.736225E-7,
     * 3.442177E-7, 2.671323E-6, 1.832858E-6,-2.394688E-7,-3.478444E-7,
     * 8.702406E-6,-8.421214E-6,-1.372341E-6,-1.455234E-6,-4.998479E-8,
     *-1.488378E-6,-1.251789E-5, 5.226868E-7,-2.049301E-7,0.E0,
     *-8.043059E-6,-2.991300E-6, 1.473654E-7,-3.154542E-7,0.E0,
     * 3.699128E-6,-3.316126E-6, 2.901257E-7, 3.407826E-7,0.e0,
     * 2.550120E-6,-1.241123E-6, 9.901116E-8, 2.210482E-7,0.e0,
     *-6.351059E-7, 2.341650E-6, 1.061492E-6, 2.878231E-7,0.e0/

	DATA CCSEC/ 1.289600E-6,5.550147E-1,2.076942E0,
     *3.102810E-5,4.035027E0, 3.525565E-1,
     *9.124190E-6,9.990265E-1,2.622706E0,
     *9.793240E-7,5.508259E0, 1.559103E1/,CCSEC3/-7.757020E-8/
     
	DATA CCKM/3.122140E-5/
	DATA DCARGM/5.1679830D0,8.3286911095275D3,
     *5.4913150D0,-7.2140632838100D3,
     *5.9598530D0, 1.5542754389685D4/
	DATA CCAMPM/
     * 1.097594E-1,2.896773E-7,5.450474E-2, 1.438491E-7,
     *-2.223581E-2,5.083103E-8,1.002548E-2,-2.291823E-8,
     * 1.148966E-2,5.658888E-8,8.249439E-3, 4.063015E-8/

C  EXECUTION

	DT=(DJE-DCT0)/DCJUL
	T=DT
	DTSQ=DT*DT
	DO 100 K=1,4
	DLOCAL=DMOD(DCFEL(1,K)+DT*DCFEL(2,K)+DTSQ*DCFEL(3,K),DC2PI)
	IF(K.EQ.1) DML=DLOCAL
100   IF(K.NE.1) F(K-1)=DLOCAL
	E=AMOD(CCSEL(1)+T*CCSEL(2)+sngl(dTSQ)*CCSEL(3),CC2PI)
	DO 300 K=1,4
300	SN(K)=SIN(AMOD(CCSEC(2,K)+T*CCSEC(3,K),CC2PI))
	PERTL=   CCSEC(1,1)*SN(1)+CCSEC(1,2)*SN(2)+
     +(CCSEC(1,3)+T*CCSEC3)*SN(3)+CCSEC(1,4)*SN(4)
	PERTR =0.
	DO 400 K=1,15
	A=DMOD(DCARGS(1,K)+DT*DCARGS(2,K),DC2PI)
	COSA=COS(A)
	SINA=SIN(A)
	PERTL=PERTL+CCAMPS(1,K)*COSA+CCAMPS(2,K)*SINA
400	PERTR=PERTR+CCAMPS(3,K)*COSA+CCAMPS(4,K)*SINA
	G=F(1)
	ESQ=E*E
	TWOE=E+E
	TWOG=G+G
	PHI=TWOE*((1.-ESQ*.125)*SIN(G)+E*.625*SIN(TWOG)+
     +     ESQ*.5416667*SIN(G+TWOG))
	DPSI=(DC1-ESQ)/(DC1+E*COS(G+PHI))

	D1PDRO=DC1+PERTR
	DML=(-(17.1822+.01742*T)*sin(f(2)-f(3))-1.3187*dsin(DML+DML)+
     * .2062*sin(2*(f(2)-f(3)))-.2274*sin(f(2)*2)+.1426*dsin(f(1))
     *-20.490)/206265.+DML
c +.0712*dsin(l)-.0517*sin(f(1)+2*DML)-.0386*sin(f(2)+f(3))-.0301*sin(l+2*f(2))
	DTL=DMOD(DML+PHI+PERTL,DC2PI)
	adlu=dmod(dtl+dc2pi*1.5,dc2pi)
	DSINLS=DSIN(DTL)
	DCOSLS=DCOS(DTL)

	PERTL=0.
	PERTP=0.
	DO 500 K=1,3
	A=DMOD(DCARGM(1,K)+DT*DCARGM(2,K),DC2PI)
	SINA=SIN(A)
	COSA=COS(A)
	PERTL=PERTL+CCAMPM(1,K)*SINA
500	PERTP=PERTP+CCAMPM(3,K)*COSA

	TL=F(2)+PERTL
	SINLM=SIN(TL)
	COSLM=COS(TL)
	SIGMA=CCKM/(1.+PERTP)
	DR=DPSI*D1PDRO
	FLATM=CCIM*SIN(F(3))
	A=SIGMA*COS(FLATM)
	DXH=DR*DCOSLS-A*COSLM
	DYH=DR*DSINLS-A*SINLM
	DZH=-SIGMA*SIN(FLATM)
	DEPS=DMOD(DCEPS(1)+DT*DCEPS(2)+DTSQ*DCEPS(3),DC2PI)
	deps=deps+(9.2025*cos(f(2)-f(3))+.5736*dcos(DML+DML)
     *+.0977*dcos(f(2)+f(2))-.0895*dcos(2*(f(2)-f(3))))/206264.8d0
c  +.0224*cos(f(1)+2*DML)+.0200*cos(f(2)+f(3))+.0129*cos(l+2*f(2))
	DCOSEP=DCOS(DEPS)
	DSINEP=DSIN(DEPS)
	sep=dsinep
	DYAH=DCOSEP*DYH-DSINEP*DZH
	DZAH=DSINEP*DYH+DCOSEP*DZH
	dr=dsqrt(dxh*dxh+dyh*dyh+dzh*dzh)
	r=dr
	adec=dasin(-dzah/dr)
	dra=dmod(datan2(-dyah,-dxh)+dc2pi,dc2pi)
	ra=dra
	eqt=dmod(DML-dra,dc2pi/2.)
	IF(ABS(EQT).GT.1.5)  EQT = EQT - SIGN(3.141593,EQT)
	GST = DMOD((dje-2451545.5D0)*dc2pi+DML,dc2pi)
	IF(GST.LT.0.) GST = GST + 6.2831853
	RETURN
	END
