!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  PROGRAM Kontur
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  (c) 2005, Christian B. Lang, Mathematische Methoden in der Physik
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  Das Programm zeichnet die Konturen gleicher Temperatur (oder gleichen
!  Potentials) fuer eine Platte der Groesse 5 x 8
!  Da das Programm keine Bildschirmpixel direkt ansteuern kann, werden
!  statt dessen Grafiksymbole gezeichnet
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

PROGRAM Kontur
!  Plattenausdehnung: x : 0 ... 3  (Pixel : 0 ... 3*90=270)
!                     y : 0 ... 5  (Pixel : 0 ... 5*90=450)}
  
   implicit none
   integer, parameter:: xs=90, ys=90          !Skalierung
   integer, parameter:: xb=3*xs, yb=5*ys      ! Pixel
   integer, parameter:: xa=0*xs, ya=0*xs      ! Pixel
  
   real, dimension(0:xb):: told,tnew
   real::tt, fun
   integer::x,y
   character (len=256)::gnucmd,gnudat
  
   gnudat='kontur_plot.dat'
   gnucmd='kontur_plot.cmd'
  
   print *,'Ich erzeuge zwei Files:'
   print *,trim(gnudat)
   print *,trim(gnucmd)
   print *,'Falls diese Files schon existieren, erfolgt eine Fehlermeldung.'
   open(unit=10,file=trim(gnudat),status='new')
   open(unit=11,file=trim(gnucmd),status='new')
   write(11,'(a)') 'set title "Temperaturverteilung"'
   write(11,'(a)') 'set xlabel "x"'
   write(11,'(a)') 'set ylabel "y"'
   write(11,'(a,i5,a,i5,a)') 'set yrange [',ya/xs,":",yb/xs,']'
   write(11,'(a,i5,a,i5,a)') 'set xrange [',xa/xs,":",xb/xs,']'
   write(11,'(a)') 'plot "'//trim(gnudat)//'" using 1:2 with points ps 0 t "Isothermen"'
   write(11,'(a)') 'pause -1'
   write(11,'(a)') 'quit'
   close(11)
   
!  Zeile 0:

   do x=xa,xb
     told(x)=int(fun(real(x)/xs,real(ya)/xs))
   end do
   
!  Restliche Zeilen:
   do y=ya+1,yb
!    Zeile berechnen:
     do x=xa,xb
        tnew(x)=int(fun(real(x)/xs, real(y)/ys))
     end do
     do x=xa+1,xb-1
!      Hat sich der Gradbereich geaendert?
       tt= tnew(x)
       if (tt .ne. tnew(x - 1)) write(10,*)real(x)/xs, real(y)/ys
       if (tt .ne. tnew(x + 1)) write(10,*)real(x+1)/xs, real(y)/ys
       if (tt .ne. told(x))     write(10,*)real(x)/xs,real(y-1)/ys
     end do
     told=tnew
   end do

   close(10)

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

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

   print *,'Berechnung und Darstellung sind abgeschlossen.'

END


REAL FUNCTION fun(x,y)
!  Temperaturverteilung
   implicit none
   real, intent(in)::x,y

   fun= 100+x*y*(y-x)
   
END FUNCTION fun



