!                                                     B Bunk 4/2001
! Tabelle von Quadrat- und Kubikwurzeln berechnen     rev    3/2012
!     und in File schreiben

module global

  real(4), parameter            :: tol = 1.e-6    ! Fehlertoleranz

end module global

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

program wurzeln

  integer, parameter                  :: nmax = 100
  real(4), dimension(nmax,3)          :: tabelle

  print*,'x0, dx, Anzahl?'            ! Daten eingeben und pruefen
  read*, x0, dx, n
  if (n > nmax) then
     print*,'Fehler: Anzahl zu gross fuer Tabelle'
     stop
  endif

  do i=1,n          ! Schleife ueber Tabellenzeilen

     x = x0 + (i-1)*dx

     tabelle(i,1) = x        ! Tabellenzeile berechnen
     tabelle(i,2) = w2(x)    ! Funktionsaufruf

     call sub3(x,w3)         ! Aufruf von Subroutine
     tabelle(i,3) = w3

     print*, tabelle(i,1:3)  ! Kontrollausdruck

  enddo             ! Ende Tabellenzeilen

  open(10,file='wurzeln.tab')         ! File oeffnen
  do i=1,n
     write(10,*) tabelle(i,1:3)       ! Zeile Schreiben
  enddo
  close(10)                           ! File schliessen
  stop
end program wurzeln

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

function w2(z)

  ! w2 <- sqrt(z) berechnen

  use global

  if (z < 0.) then                    ! Argument pruefen
     print*,'Fehler in w2: negatives Argument'
     stop
  endif

  w2 = 1.                                ! Startwert
  do
     w2alt = w2
     w2 = .5 * (w2 + z/w2)               ! Iteration
     if ( abs(w2-w2alt) <= tol ) exit    ! Konvergenztest
  enddo
end function w2

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine sub3(z,w3)

  ! w3 <- z**(1/3) berechnen

  use global

  if (z < 0.) then                    ! Argument pruefen
     print*,'Fehler in w3: negatives Argument'
     stop
  endif

  w3 = 1.                                ! Startwert
  do
     w3alt = w3
     w3 = (2*w3 + z/w3**2)/3             ! Iteration
     if ( abs(w3-w3alt) <= tol ) exit    ! Konvergenztest
  enddo
end subroutine sub3