! Berechne die Flaeche, die von Enden der Zeiger einer Analoguhr aufgespannt wird

program analoguhr
  real,parameter::pi=4*atan(1.0)
  real,dimension(3)::lambda
  real::t  
  integer::i,n,nn,nmax

  !lambda=(/0.5,1.0,0.8/)
  print *,"Lambda 1:"; read *,lambda(1)
  print *,"Lambda 2:"; read *,lambda(2)
  print *,"Lambda 3:"; read *,lambda(3)

  ! Suche lokale Maxima
  call flaeche_max_int(lambda,6*3600)
  call flaeche_max(lambda)
  stop

contains
  function uhrzeit(t)
      real::t
      character(len=17)::uhrzeit,str
    integer::h,m,i
    real::s
    h=floor(t)
    m=floor((t-h)*60)
    s=((t-h)*60-m)*60
    write (unit=str,fmt="(a,i2,a,i2,a,f9.6,a)") "[",h,":",m,":",s,"]"
    do i=1,len(str)
      if (str(i:i)==" ") str=str(:i-1)//"0"//str(i+1:)
    enddo
    uhrzeit=str
  end function
  function flaeche(lambda,t)          ! Berechnet die Flaeche zur Zeit t (in Stunden)
      real,dimension(3)::lambda
      real::t,flaeche
    real::phi1,phi2,phi3,s12,s23,s31
    phi1=pi*t/6.0
    phi2=pi*t/6.0*12
    phi3=pi*t/6.0*12*60
    s12=lambda(1)*lambda(2)*sin(phi2-phi1)
    s23=lambda(2)*lambda(3)*sin(phi3-phi2)
    s31=lambda(3)*lambda(1)*sin(phi1-phi3)
    flaeche=0.5*(s12+s23+s31)
  end function
  function flaeche_t(lambda,t)         ! Ableitung der Flaeche zur Zeit t (in Stunden)
      real,dimension(3)::lambda
      real::t,flaeche_t
    real::phi1,phi2,phi3,s12,s23,s31
    phi1=pi*t/6.0
    phi2=pi*t/6.0*12
    phi3=pi*t/6.0*12*60
    s12=lambda(1)*lambda(2)*(phi2-phi1)/t*cos(phi2-phi1)
    s23=lambda(2)*lambda(3)*(phi3-phi2)/t*cos(phi3-phi2)
    s31=lambda(3)*lambda(1)*(phi1-phi3)/t*cos(phi1-phi3)
    flaeche_t=0.5*(s12+s23+s31)
  end function
  subroutine flaeche_max_int(lambda,n) ! Berechnet die maximalen Flaechen, n=Teilintervalle
      real,dimension(3)::lambda
      integer::n
    integer::i,j
    real::t,t1,t2,f,f_t1,f_t2,t_max,f_max
    t1=0.0; f_t1=flaeche_t(lambda,t1)
    t_max=0.0; f_max=0.0
    do i=1,n
      t2=6.0*i/real(n); f_t2=flaeche_t(lambda,t2)
      if (f_t1*f_t2<0) then            ! Vorzeichenwechsel der Ableitung=Maximum der Fleache
        do j=1,50                      ! Bisektionssuche
          t=0.5*(t1+t2)
          if (flaeche_t(lambda,t)*f_t1>0) then
            t1=t
          else 
            t2=t
          endif
        end do
        f=flaeche(lambda,t)
        if (abs(f)>f_max) then
          f_max=abs(f); t_max=t
        end if
        print *,"Lokales Maximum f=",f, "  bei t=",t/6.0,uhrzeit(t)
      end if
      t1=t2; f_t1=f_t2
    end do
    print *,"Globales Maximum f=",f_max, "  bei t=",t_max/6.0,uhrzeit(t_max)
  end subroutine
  subroutine flaeche_max(lambda)       ! Berechnet die maximalen Flaechen analytisch
      real,dimension(3)::lambda
    real::a,b,d,x,p,q,dis,r,phi,f
    p=-sum(lambda**2)/(lambda(1)*lambda(2))**2  ! Loese x^3+px+q=0
    q=-2.0*lambda(3)/(lambda(1)*lambda(2))**2
    dis=(q/2.0)**2+(p/3.0)**3          ! Cardano-Formeln
    r=sqrt(-dis+(q/2.0)**2)            ! (r,phi)=Polarkoordinaten von u^3
    phi=atan(sqrt(-dis)/(-q/2.0))
    x=2.0*(r**(1.0/3.0))*cos(phi/3.0)  ! x=u+v
    d=1.0/x
    a=sqrt(lambda(1)**2-d**2)  
    b=sqrt(lambda(2)**2-d**2)
    f=0.5*(a+b)*(d+lambda(3))
    print *,"Maximale Flaeche (freie Zeiger):",f
  end subroutine

end program