Gittata pale eoliche

'-------------------------------------------------------
' Distanza di caduta e tempo di volo di un corpo
' lanciato in aria, in presenza di attrito viscoso.
'-------------------------------------------------------

' Variabili globali
Dim V_x As Double, V_y As Double ' componenti velocità
Dim Gamma As Double ' coefficiente attrito
Dim Acc_g As Double ' acc. gravità in m/s^2

' Costanti
Const GR_RD As Double = 1.74532925199433E-02 ' gradi -> rad


Public Sub Calcola_Gittata()
   Dim v_0 As Double ' velocita iniziale
   Dim alpha As Double ' angolo di lancio
   Dim t_star As Double
   
   ' Accelerazione di gravità
   Acc_g = Range("B36")

   ' Coefficiente di attrito
   Gamma = Range("L40")
   If Gamma < 0# Then
       MsgBox "Attrito negativo!", vbExclamation & vbOKOnly, "Errore!"
       Exit Sub
   End If
   If Gamma < 0.000001 Then
       Gamma = 0.000001
       Range("L40") = Gamma
   End If

   ' Modulo velocità
   v_0 = Range("B40")
   If v_0 <= 0# Then
       MsgBox "Il modulo velocità deve essere positivo!", vbExclamation & vbOKOnly, "Errore!"
       Exit Sub
   End If
   
   '-------------------------
   
   'Angolo di lancio
   alpha = Range("B44")
   alpha = alpha * GR_RD ' Gradi -> Radianti
   
   ' Calcolo le componenti velocità
   V_x = v_0 * Cos(alpha)
   V_y = v_0 * Sin(alpha)
   If V_y <= 0# Then
       MsgBox "Lancio verso il basso!", vbExclamation & vbOKOnly, "Errore!"
       Exit Sub
   End If
   
   ' altezza punto di partenza
   Dim H As Double
   H = Range("G44")

   ' trova t*
   t_star = tempo_impatto(H)
   Range("M44") = trova_x(t_star, -Range("E44"))
   Range("L44") = t_star
   
   '-------------------------

   'Angolo di lancio
   alpha = Range("B48")
   alpha = alpha * GR_RD ' Gradi -> Radianti
   
   ' Calcolo le componenti velocità
   V_x = v_0 * Cos(alpha)
   V_y = v_0 * Sin(alpha)
   If V_y <= 0# Then
       MsgBox "Lancio verso il basso!", vbExclamation & vbOKOnly, "Errore!"
       Exit Sub
   End If
   
   ' altezza punto di partenza
   H = Range("G48")
       
   ' trova t*
   t_star = tempo_impatto(H)
   Range("M48") = trova_x(t_star, Range("E48"))
   Range("L48") = t_star
   
End Sub

'-------------
' tempo_impatto()
'
' tempo di impatto con il metodo della secante
' PRECISIONE : Precisione richiesta
'-------------
Private Function tempo_impatto(H As Double) As Double

Const PRECISIONE As Double = 0.00000001
Const N_MAX_ITERAZIONI As Integer = 1000
Dim t1 As Double
Dim t2 As Double
Dim t3 As Double
Dim Y_t1 As Double
Dim Y_t2 As Double
Dim Y_t3 As Double
Dim i As Integer
Dim denom As Double

   ' assegno valori iniziali
   t1 = 0#
   t2 = 100#
   Y_t1 = trova_y(t1, H)
   Y_t2 = trova_y(t2, H)
   
   For i = 1 To N_MAX_ITERAZIONI
       denom = (Y_t2 - Y_t1)
       If denom = 0# Then Exit For
       t3 = (t1 * Y_t2 - t2 * Y_t1) / denom
       Y_t3 = trova_y(t3, H)
       If Abs(Y_t3) < PRECISIONE Then
           tempo_impatto = t3
           Exit Function
       End If
       If Sgn(Y_t3) <> Sgn(Y_t1) Then
           t2 = t3
           Y_t2 = Y_t3
       Else
           t1 = t3
           Y_t1 = Y_t3
       End If
   Next i
   
   ' Precisione troppo alta
   tempo_impatto = 0#

End Function

Private Function trova_y(t As Double, H As Double) As Double
   trova_y = -(Acc_g / Gamma) * t + (1# / Gamma ^ 2) * Acc_g + (1# / Gamma) * V_y + exp(-Gamma * t) * (H - (V_y / Gamma) - (Acc_g / (Gamma ^ 2)))
End Function

Private Function trova_x(t As Double, X0 As Double) As Double
   trova_x = (V_x / Gamma) + (X0 - V_x / Gamma) * exp(-Gamma * t)
End Function

'-------------
Public Sub Trova_angoli()
   Dim alfa As Double
   alfa = angolo_gittata_max("B44", "J44")
   Range("B44") = alfa
   alfa = angolo_gittata_max("B48", "J48")
   Range("B48") = alfa
End Sub

'-------------
' angolo_gittata_max()
'
' PRECISIONE : Precisione richiesta
'-------------
Private Function angolo_gittata_max(S1 As String, S2 As String) As Double

Const PRECISIONE As Double = 0.0001
Const N_MAX_ITERAZIONI As Integer = 100
Dim alfa1 As Double
Dim alfa2 As Double
Dim g_alfa1 As Double
Dim g_alfa2 As Double
Dim i As Integer

   ' assegno valori iniziali
   alfa1 = 0#
   alfa2 = 90#
   g_alfa1 = trova_g(alfa1, S1, S2)
   g_alfa2 = trova_g(alfa2, S1, S2)

   For i = 1 To N_MAX_ITERAZIONI
       If g_alfa1 > g_alfa2 Then
           alfa2 = alfa2 + (alfa1 - alfa2) / 2
           g_alfa2 = trova_g(alfa2, S1, S2)
       Else
           alfa1 = alfa1 + (alfa2 - alfa1) / 2
           g_alfa1 = trova_g(alfa1, S1, S2)
       End If
   
       If Abs(alfa1 - alfa2) < PRECISIONE Then
           angolo_gittata_max = (alfa1 + alfa2) / 2
           Exit Function
       End If
   Next i
   
   ' Precisione troppo alta
   angolo_gittata_max = 45#

End Function

Private Function trova_g(alfa As Double, S1 As String, S2 As String) As Double
   Range(S1) = alfa
   trova_g = Range(S2)
End Function