Program pohlschesrad
!	use dislin
	Implicit None

	!Definitionen der Variablen
	Integer, parameter :: tmax = 3100	!Dargestelltes Zeitintervall
	Integer, parameter :: tmin = 3020
	Real(8), parameter :: dt = 0.001		!Zeitschritt
	Integer, parameter :: n  = int(real(tmax)/dt)	!Anzahl der Schritte bei der Integration
	Integer, parameter :: m  = int(real(tmin)/dt)	!Die ersten m Schritte werden nicht dargestellt (Einschwingverhalten)
	Integer, parameter :: mm = 100		!Die ersten mm Minima werden in der Return-Map nicht dargestellt ( " )
	

	Real(8) :: x(n)				!Phi
	Real(8) :: y(n)				!Phipunkt
	Real    :: it(500)			!Array der Minima fuer die Return-Map
	Real(8) :: k1(2), k2(2), k3(2), k4(2)	!Runge-Kutta-Parameter, 1. Komp. x, 2. Komp. y
	Real(8) :: dpf, nl, antr, rsk, om	!Systemparameter

	!Variabln zur graphischen Darstellung
	Real(8) :: zeit(n-m)
	Real(8) :: xpl(n-m)
	Real(8) :: ypl(n-m)
	Real, allocatable :: itpl1(:), itpl2(:)
	Real :: d1,d2,d3,d4,d5,d6,d7,d8
	Integer :: mmm

	Integer :: i, j				!Indizes fuer Schleifen

	Real(8), parameter :: sechstel = 1.0d0/6.0d0 !Spart Rechenzeit


	!Festlegen der Systemparameter

	rsk 		= 9.44		!Direktionskonstante
	dpf		= 0.799		!Daempfungskonstante
	nl 		= 14.68		!Nichtlinearitaet
	antr 		= 2.1		!Antriebsamplitude
	om		= 2.178 		!Kreisfrequenz des Antriebs


	!Festlegen der Anfangswerte

	x(1)=0
	y(1)=0


	!Klassisches Runge Kutta-Verfahren 4.Ordnung

	Do i=1, n-1

		k1(1)= f1( x(i) , y(i) ) * dt
		k1(2)= f2( x(i) , y(i) , t(i) ) * dt

		k2(1)= f1( x(i) + 0.5d0*k1(1) , y(i) + 0.5d0*k1(2) ) * dt
		k2(2)= f2( x(i) + 0.5d0*k1(1) , y(i) + 0.5d0*k1(2) , t(i) + 0.5d0*dt )* dt

		k3(1)= f1( x(i) + 0.5d0*k2(1) , y(i) + 0.5d0*k2(2) ) * dt
		k3(2)= f2( x(i) + 0.5d0*k2(1) , y(i) + 0.5d0*k2(2) , t(i) + 0.5d0*dt) * dt

		k4(1)= f1( x(i) + k3(1), y(i) + k3(2) ) * dt
		k4(2)= f2( x(i) + k3(1), y(i) + k3(2), t(i) + dt ) * dt

		x(i+1)= x(i) + sechstel * ( k1(1) + 2.0d0 * k2(1) + 2.0d0 * k3(1) + k4(1))
		y(i+1)= y(i) + sechstel * ( k1(2) + 2.0d0 * k2(2) + 2.0d0 * k3(2) + k4(2))

	End Do


	!Berechnung der Return-Map

	j=2

	Do i=1,n-1

		If ( y(i) <= 0 .And. y(i+1) >= 0 ) Then
			it(j) = real(x(i))
			j=j+1
		End If

	End Do



	mmm = j - 1 - mm
	allocate(itpl1(j-2))
	allocate(itpl2(j-2))

	do i=mm,j-2
		itpl1(i) = 3.1414/2. - it(i)
		itpl2(i) = 3.1414/2. - it(i+1)
	end do
	do i=m+1,n
		zeit(i-m) = i*dt
		xpl(i-m) = x(i)
		ypl(i-m) = y(i)
	end do


	call scrmod('revers')
	!call metafl ('png')
	call filmod('delete')

	!call metafl('png')
	call disini ()

	call texmod('on')

	call messag('$\partial_t^2 \phi + \frac{K}{J} \, \, \partial_t \phi + \frac{D}{J} \, \phi - \frac{N}{J} \, \sin \phi = \frac{F}{J} \, \sin(\omega t$)',650,100)
	call messag('$\frac{K}{J} =$',2400,200)
	call messag('$\frac{D}{J} =$',2400,350)
	call messag('$\frac{N}{J} =$',2400,500)
	call messag('$\frac{F}{J} =$',2400,650)
	call messag('$\Large{\omega} =$',2400,800)
	call number(Real(dpf),3,2475,200)
	call number(Real(rsk),3,2475,350)
	call number(Real(nl),3,2475,500)
	call number(Real(antr),3,2475,650)
	call number(Real(om),5,2500,800)

	call titlin ('Return-Map',4)
	call axspos(1450,900)
	call axslen(800,600)
	!call setscl(real(itpl1),mm,'x')
	!call setscl(real(itpl2),mm,'y')
	call name ('$\pi/2-\phi_{min,N}','x')
	call name ('$\pi/2-\phi_{min,N+1}','y')
	call graf(0.8,1.6,0.8,.4,0.8,1.6,0.8,.4)
	call title
	call incmrk(-1)
	call hsymbl(1)
	call curve (itpl1,itpl2,mmm)
	call endgrf


	!call gaxpar(0.,3.,'extend','x',xmin,xmax,ixor,xstp,xndig)
	!call gaxpar(-3.,3.,'extend','y',ymin,ymax,yor,ystp,yndig)

	call incmrk(0)
	call titlin ('Phasenraum',4)
	call axspos(300,900)
	call axslen(800,600)
	call setscl(real(xpl),n-m,'x')
	call setscl(real(ypl),n-m,'y')
	call name ('$\phi$','x')
	call name ('$\frac{d\phi}{dt}$', 'y')
	call graf(d1,d2,d3,d4,d5,d6,d7,d8)
	call title
	call curve(real(xpl),real(ypl),n-m)
	call endgrf


	call titlin ('Zeitverlauf', 4)
	call setscl(real(zeit),n-m,'x')
	call setscl(real(xpl),n-m,'y')
	call axspos(250,1800)
	call axslen(2500,600)
	call name ('$t$','x')
	call name ('$\phi(t)$','y')
	call graf(d1,d2,d5,d6,d3,d4,d7,d8)
	call title
	call curve(real(zeit),real(xpl),n-m)
	call endgrf

	call disfin ()
	call symfil ('cons','delete')

!Definition der Funktionen

contains

Real(8) Function t(i)
	Implicit None
	integer i
	t=i*dt
End Function

Real(8) Function f1(x,y)
	Implicit None
	Real*8 x, y
	f1=y
End Function

Real(8) Function f2(x,y, t)
	Implicit None
	Real*8 x, y, t
	f2= -rsk*x -dpf*y +antr*Cos(om*t) + nl*sin(x)
End Function


End Program