!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  PROGRAM Rotate
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  (c) 2005, Christian B. Lang, Mathematische Methoden in der Physik
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  Animation eines Drahtgittermodells. Das Programm stellt eine gegebene
!  Menge von Strecken dreidimensional dar. Weiters wir das gegebene Objekt
!  im Raumgedreht und so "animiert". Der Blickwinkel des Betrachters und
!  die Drehwinkel werden erfragt. Die Punkte und Strecken sind in einem
!  vorgespeichert und sollten bei Bedarf geaendert oder eingelesen werden.
!  Vorgeschlagene Testeingaben:
!  phi, theta: 1., 1.
!  dphi. dtheta, dpsi:  -1.57, 0.1 ,1.57 (um y-Achse)
!  oder                    .1, 0   , 0   (um z-Achse)
!  oder                    0 , 0.1 , 0   (um x-Achse)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

PROGRAM Rotate
   implicit none
   
!  Testdaten: Punkte
   integer, parameter:: np = 8
   real, dimension(1:3,1:np)::ptest
   parameter ( ptest = &
   &  (/(/1.,1.,-1./), (/-1.,1.,-1./), (/-1.,-1.,-1./), (/1.,-1.,-1./), &
   &    (/1.,1.,1./), (/-1.,1.,1./), (/-1.,-1.,1./), (/1.,-1.,1./)/) )

!  Testdaten: Strecken
   integer, parameter:: ns = 12
   integer, dimension(1:2,1:ns)::stest
   parameter ( stest = &
   &  (/ (/1,2/), (/2,3/), (/3,4/),  (/4,1/),&
   &     (/5,6/), (/6,7/), (/7,8/),  (/8,5/),&
   &     (/1,5/), (/2,6/), (/3,7/),  (/4,8/)/))

   real::csa,csb,csc,sna,snb,snc,dphi,dtheta,dpsi
   real, dimension(1:3,1:3)::r
   integer, dimension(1:2,1:ns)::s
   real, dimension(1:3,1:np)::ps,p
   integer::i
   character (len=1)::antwort
   character (len=256)::gnucmd,gnuloop,gnudat

   s=stest
   p=ptest
  
   write(*,*) 'Drehung pro Zeitschritt um dphi, dtheta, dpsi (in rad):'
   write(*,*) 'dphi   (z.B. 0.05):'
   read(*,*) dphi
   write(*,*) 'dtheta (z.B. 0.05):'
   read(*,*) dtheta
   write(*,*) 'dpsi   (z.B. 0   ):'
   read(*,*) dpsi

!  Drehmatrix berechnen
   csa = COS(dphi)
   sna = SIN(dphi)
   csb = COS(dtheta)
   snb = SIN(dtheta)
   csc = COS(dpsi)
   snc = SIN(dpsi)
   r(1, 1) = csc * csa - csb * sna * snc
   r(1, 2) = csc * sna + csb * csa * snc
   r(1, 3) = snb * snc
   r(2, 1) = -csa * snc - csb * sna * csc
   r(2, 2) = -sna * snc + csb * csa * csc
   r(2, 3) = snb * csc
   r(3, 1) = snb * sna
   r(3, 2) = -snb * csa
   r(3, 3) = csb
   print *,r

   gnudat='rotate_plot.dat'
   gnuloop='rotate_loop.cmd'
   gnucmd='rotate_plot.cmd'

   print *,'Ich erzeuge drei Files:'
   print *,trim(gnudat)
   print *,trim(gnuloop)
   print *,trim(gnucmd)
   print *,'Falls diese Files schon existieren, erfolgt eine Fehlermeldung.'
   open(unit=9 ,file=trim(gnuloop),status='new')
   open(unit=10,file=trim(gnudat),status='new')
   open(unit=11,file=trim(gnucmd),status='new')

   write(9,'(a)') 'plot "'//trim(gnudat)//'" with lines'
   write(9,'(a)') 'pause 0.2'
   write(9,'(a)') 'if(time>0)reread'
   close(9)

   write(11,'(a)') 'set title "Rotate"'
   write(11,'(a)') 'set xrange [-2:2]'
   write(11,'(a)') 'set yrange [-2:2]'
   write(11,'(a)') 'set size ratio 1'
   write(11,'(a)') 'time=10'
   write(11,'(a)') 'load "'//trim(gnuloop)//'" '
   write(11,'(a)') 'quit'
   close(11)

   write(*,*) 'In der endlosen Schleife wird nun die Drehung animiert.'
   write(*,*) 'Nach jedem Bild muessen Sie CR eingeben, um das naechste Bild zu erzeugen.'
   write(*,*) 'Dies kann nur durch einen Abbruch (unter Unix: CTRL-C) gestoppt werden.'
   write(*,*) 'Soll ich die Animation starten? (j oder n):'
   read(*,'(a1)')antwort

   if((antwort=='j').or.(antwort=='J'))then
     call Plot(s,ns,p)
     close(10)

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

     call system(trim('gnuplot '//trim(gnucmd)//' &'))
    
     do
       open(unit=10,file=trim(gnudat),status='unknown')
       call Plot(s,ns,p)
       close(10)
       read(*,*)
       call Drehung(p, r, ps, np)
       p=ps
!       Werte rueckspeichern
     end do
   else
     write(*,*) 'Wunschgemaess abgebrochen.'
   end if
    
END PROGRAM Rotate

SUBROUTINE Drehung (p,r,ps,npt)
!  Drehung PS = R P
   implicit none

   integer, intent(in)::npt
   real, dimension(1:3,1:npt), intent(in)::p
   real, dimension(1:3,1:npt), intent(out)::ps
   real, dimension(1:3,1:3), intent(in)::r
   integer::i

   do i=1,npt
     ps(:,i)=matmul(r(:,:),p(:,i))
   end do

END SUBROUTINE Drehung

SUBROUTINE Plot (s,nst,p)
!  Zeichnen aller Teilstrecken
   implicit none

   integer, intent(in)::nst
   integer, dimension(1:2,1:nst)::s
   real, dimension(1:3,1:1), intent(in)::p
   integer::i

   do i = 1,nst
     write(10,*)p(1, s(1, i)), p(2, s(1, i))
     write(10,*)p(1, s(2, i)), p(2, s(2, i))
     write(10,*)
   end do

END SUBROUTINE Plot


