      program krzywe_stozk
c----------------------------------------------------------------------c
c Wylicza wspolrzedne punktow do wykresow krzywych stozkowych.
c UWAGA!
c Program nalezy skompilowac przy uzyciu dowolnego kompilatora 
c Fortranu 90, np. gfortran, wykonując komendę
c     gfortran krzywe_stozk.f
c a następnie wykonac komendę
c     ./a.out
c i wykreslic wygenerowany zbior punktow o nazwie "ks.dat" przy uzyciu
c dowolnego programu do robienia wykresów, np. gnuplot.
c Rozne krzywe stozkowe otrzymamy dla roznych wartosci mimosrodu "eps".
c
c Karol Kolodziej, University of Silesia, karol.kolodziej@us.edu.pl
c 26 January, 2012; modified 18 Nov., 2012, 1 May 2017
c----------------------------------------------------------------------c

      implicit none

      integer, parameter :: n=100000 ! liczba punktow
      real(kind(1.d0)), parameter :: dmax=1.d0 ! maksymalna odleglosc 
                                                 ! dwoch kolejnych x-ow
      real(kind(1.d0)) :: p,eps,twopi,phi,dphi,x,y,x0

      p=1.0d0    ! dowolny parametr
c Mimosrod:
      eps=0.0000001d0     ! okrag (prawie), eps=0
      eps=0.5d0           ! elipsa,  0 < eps < 1
      eps=1.0d0           ! parabola, eps=1
      eps=1.0002d0        ! hiperbola, eps > 1
c      eps=1.02d0        ! hiperbola, eps > 1

      twopi=8.d0*atan(1.d0)

      dphi=twopi/n

      open(unit=11,file='ks.dat',status='replace',action='write',
     &                                                       err=200)
      write(11,'(a,e15.8,a,e15.8)') '# Krzywa stozkowa z p=',p,
     &                ' i e=',eps

      phi=0.d0
      x0=r(phi)*cos(phi)

      do
       phi=phi+dphi
       if (phi > twopi) exit
       x=r(phi)*cos(phi)
       y=r(phi)*sin(phi)
       if (abs(x-x0) < dmax) then
        write(11,'(2(e16.8))') x,y
        x0=x
       end if
      end do

      close(11)

      stop

 200  write(*,*)'Something is wrog with the files'

      contains

      function r(phi)
       real(kind(1.d0)) :: r,phi
       r=p/(1.d0+eps*cos(phi))
      end function r

      end
