Fortran: Anhang B: Kettenlinie: Catmod1.f90

module catmod1
  implicit none 

  real, parameter :: PI = 3.141592654
  real            :: length, width, height, q
  real            :: a, x1, x2, y1, y2, s1, s2, phi1, phi2, V1, V2, H, F1, F2, &
                     Fmax, sag, xi1, xi
		 
  private :: seilparameter

  contains
    subroutine input
      write (*,"(/ A /)") "********** KETTENLINIE **********"
      write (*, "(A)") " Eingabe:"
      write (*, "(A, $)") "  Seillänge                  L: "
      read (*, *) length
      write (*, "(A, $)") "  Breite                     b: "
      read (*, *) width
      write (*, "(A, $)") "  Höhe                       h: "
      read (*, *) height		
      write (*, "(A, $)") "  Spezifisches Kettengewicht q: "
      read (*, *) q				
    end subroutine input
	

    ! Bestimmung des Seilparameters u.a. mittels Newton-Näherungsverfahren
    subroutine seilparameter()
      real               :: func, func_der, k, xi_old=-999.9
      integer            :: i = 0
      integer, parameter :: i_max = 100

      if (width<=epsilon(width)) then
        write(*,*) "Ungültiger Wert für die Breite b! Programmabbruch!"
        stop
      end if 

      if ((length**2-width**2-height**2) <= epsilon(length)) then
        write(*,*) "Seillänge L zu kurz! Programmabbruch!"
        stop			
      end if			
	
      k = sqrt(length**2 - height**2) / width
  		
      ! 1. Näherung für Startwert (aus Reihenentwicklung)
      xi1 = sqrt(-10.0 + sqrt(120.0*k-20.0))  
      xi = xi1
	
      ! Newton-Näherungsverfahren
      do while (abs(xi_old-xi) >= epsilon(xi) .AND. i<=i_max)
        i = i+1
        xi_old = xi

        func = k*xi_old - sinh(xi_old)
        func_der = k - cosh(xi_old)

        xi = xi_old - func/func_der
      end do		

      a = width / (2*xi)
    end subroutine seilparameter


    subroutine calculate()
      real :: k
			
      call seilparameter
		
      s1 = height/2.0 * (exp(2*xi)+1)/(exp(2*xi)-1) - length/2.0
      s2 = s1 + length
      phi1 = atan2(s1, a)
      phi2 = atan2(s2, a)
      x1 = a * log(s1/a + sqrt((s1/a)**2+1))
      x2 = x1 + width
      y1 = sqrt(a**2 + s1**2)
      y2 = y1 + height
      V1 = q * s1
      V2 = q * s2
      F1 = q * sqrt(a**2 + s1**2)
      F2 = q * sqrt(a**2 + s2**2)
      H = sqrt(F1**2 - V1**2)
		
      k = height/width
      sag = y1 + k*(a*log(k + sqrt((k)**2+1))-x1) - a*sqrt(k**2+1)
    end subroutine calculate
	
	
    subroutine output
      write (*, "(/ A )") "*** Ergebnisse ***"
      write (*,*) " Eingabe:"
      write (*,*) "  L                            = ", length
      write (*,*) "  b                            = ", width
      write (*,*) "  h                            = ", height
      write (*,*) "  q                            = ", q
      write (*,*) " Berechnet:"                 
      write (*,*) "  Seilparameter a (1.Näherung) = ", width/(2*xi1)
      write (*,*) "  Seilparameter a              = ", a
      write (*,*) "  x1                           = ", x1
      write (*,*) "  x2                           = ", x2
      write (*,*) "  y1                           = ", y1
      write (*,*) "  y2                           = ", y2
      write (*,*) "  s1                           = ", s1
      write (*,*) "  s2                           = ", s2
      write (*,*) "  phi1                         = ", phi1, "rad"
      write (*,*) "  phi1                         = ", phi1 * 180.0 /PI, "°"
      write (*,*) "  phi2                         = ", phi2, "rad"
      write (*,*) "  phi2                         = ", phi2 * 180.0 /PI, "°"
      write (*,*) "  F1                           = ", F1
      write (*,*) "  F2                           = ", F2
      write (*,*) "  H=H1=H2                      = ", H
      write (*,*) "  V1                           = ", V1
      write (*,*) "  V2                           = ", V2
      write (*,*) "  Durchhang f                  = ", sag
    end subroutine output
end module catmod1