!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  PROGRAM PartSum
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  (c) 2005, Christian B. Lang, Mathematische Methoden in der Physik
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  Die Partialsummen S_n(x) (festes n) werden fuer
!  verschiedene x-Werte mit der durch die Potenzreihe
!  dargestellten Funktion f(x) verglichen. Der Wert der
!  Partialsumme wird in der Funktion S_n(x) berechnet. }
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

PROGRAM PartSum

   implicit none

   integer:: n,i
   real:: a,b,xa,xb,x,x_del,x_old,f,Reihe
   character (len=256)::gnucmd,gnudat1,gnudat2

   write(*,*) 'Ich berechne die Partialsumme S_n(x) und vergleiche ',&
   &        'mit der Funktion f(x).'
   write(*,*) 'Bitte geben Sie einen Wert fuer n an {z.B.: 5}: '
   read(*,*) n
   write(*,*) 'Ich berechne die Werte zwischen xa und xb.'
   write(*,*) 'xa = {z.B.: -9}'
   read(*,*) xa
   write(*,*) 'xb = {z.B.: 9}'
   read(*,*) xb
   write(*,*) 'Schranken fuer die Darstellung der Summe und Funktionswerte:'
   write(*,*) 'Untere Schranke = {z.B.:-4}'
   read(*,*) a
   write(*,*) 'Obere  Schranke = {z.B.: 4}'
   read(*,*) b

   gnudat1='partsum_plot1.dat'
   gnudat2='partsum_plot2.dat'
   gnucmd='partsum_plot.cmd'

   print *,'Ich erzeuge drei Files:'
   print *,trim(gnudat1)
   print *,trim(gnudat2)
   print *,trim(gnucmd)
   print *,'Falls diese Files schon existieren, erfolgt eine Fehlermeldung.'
   open(unit=9 ,file=trim(gnudat1),status='new')
   open(unit=10,file=trim(gnudat2),status='new')
   open(unit=11,file=trim(gnucmd),status='new')
   write(11,'(a)') 'set title "Partialsummen"'
   write(11,'(a)') 'set xlabel "x"'
   write(11,'(a)') 'set ylabel "y"'
   write(11,'(a,f5.1,a,f5.1,a)') 'set yrange [',a,":",b,']'
   write(11,'(a)') 'plot "'//trim(gnudat1)//'" using 1:2 with linespoints pt 4 t "Reihe", "'&
   &                      //trim(gnudat2)//'" using 1:2 with linespoints pt 2 t "Funktion"'
   write(11,'(a)') 'pause -1'
   write(11,'(a)') 'quit'
   close(11)

   x_del=(xb-xa)/100

!  Zuerst wird die Reihe geplottet

   write(9,*)xa,Reihe(xa,n)
   x_old=xa

   do i=1,100
      x= x_old + x_del
      write(9,*)x,Reihe(x,n)
      x_old=x
   end do

!  Nun wird die Funktion geplottet

   write(10,*)xa,f(xa)
   x_old=xa

   do i=1,100
      x= x_old + x_del
      write(10,*)x,f(x)
      x_old=x
   end do

   close(9)
   close(10)

!  Wenn das Unterprogramm "system" nicht verfuegbar ist, so kommentieren
!  Sie die folgende Zeile aus. Dann muessen Sie die Grafik selbst mittels
!  "gnuplot partsum_plot.cmd" aufreufen.

   call system(trim('gnuplot '//trim(gnucmd)))

   print *,'Berechnung und Darstellung sind abgeschlossen.'

END PROGRAM partsum

REAL FUNCTION f(x)

   implicit none
   real :: x
  
   f = cos(x)

END FUNCTION f

REAL FUNCTION Reihe(x,n)

   implicit none
   real, intent(in) :: x
   integer, intent(in):: n
   real :: ck,Summe
   integer:: k

   Summe = 1.0

   if(n > 1)then
     ck = 1
     do k = 1, n
       ck = - ck*x*x / (2*k*(2*k-1))
       Summe = Summe + ck
     end do
   end if
  
   Reihe=Summe
   
END FUNCTION Reihe

