Find an angle to launch the projectile at to reach a specific point

Question

So, Dani in his slightly new video -> "Making a Game, But I Only Have 3 Days" (https://youtu.be/S7Dl6ATRK2M) made a enemy which has a bow and arrow (at 5:39). I tried to recreate that but had no luck... I also can't find the website that he used... Today I found this https://physics.stackexchange.com/questions/56265/how-to-get-the-angle-needed-for-a-projectile-to-pass-through-a-given-point-for-t. It worked very well but still had problems if the target was far away and also it wasn't as accurate. The code so far is

float CalculateAngle()
    {
        float gravity = Physics.gravity.magnitude;
        float deltaX = targetPositionMod.x - currentPosition.x;
        float deltaY = targetPositionMod.y - currentPosition.y;
        float RHSFirstPart = (velocity * velocity) / (gravity * deltaX);
        float RHSSecondPart = Mathf.Sqrt(
                ((velocity * velocity) * ((velocity * velocity) - (2 * gravity * deltaY)) 
                                / (gravity * gravity * deltaX * deltaX))
                                                - 1);
        float tanθ = RHSFirstPart - RHSSecondPart;
        float angle = Mathf.Atan2(tanθ, 1) * Mathf.Rad2Deg;
        if (angle < 0) return angle;
        return -angle;
    }

The -angle is because the forward axis starts points up when the x-rotation is negative (Unity). Maybe the reason of this not working as intended is that I am not that good at this kind of Physics (Part of that is me being only 14). Maybe the problem is in the code, maybe it is the formula. Any help is appreciated.

Thanks...

Edit: The Archer class is:

using UnityEngine;
using System;
public class Archer : MonoBehaviour
{
    [SerializeField] float velocity = default;
    [SerializeField] Transform target = default;
    [SerializeField] GameObject arrowPrefab = default;
    [SerializeField] float coolDown = default;
    Vector3 targetPositionMod;
    Vector3 currentPosition;
    Vector3 targetPosition;
    float countDown = 0f;

    void Start()
    {
        countDown = coolDown;
        UpdateVariables();
    }

    void Update()
    {
        UpdateVariables();
        SetAngle();
        ShootBullet();
    }
    
    void UpdateVariables()
    {
        currentPosition = transform.position;
        targetPositionMod = Mod(target.position);
        targetPosition = target.position;
        targetPosition.x /= 10;
        targetPosition.y /= 10;
        targetPosition.z /= 10;
        countDown -= Time.deltaTime;
    }

    void SetAngle()
    {
        Vector3 direction = targetPosition - currentPosition;
        Quaternion lookRotation = Quaternion.LookRotation(direction);
        Vector3 rotation = lookRotation.eulerAngles;
        rotation.x = (float) CalculateAngle();
        transform.rotation = Quaternion.Euler(rotation.x, rotation.y, 0f);
    }

    void ShootBullet()
    {
        if (!(countDown <= 0f)) return;
        countDown = coolDown;
        GameObject arrow = Instantiate(arrowPrefab, transform.position, transform.rotation);
        Rigidbody Rigidbody = arrow.GetComponent<Rigidbody>();
        Rigidbody.AddForce(transform.forward * velocity, ForceMode.Impulse);
    }

    double CalculateAngle()
    {
        double gravity = Physics.gravity.magnitude;
        double deltaX = targetPositionMod.x - currentPosition.x;
        double deltaY = targetPositionMod.y - currentPosition.y;
        double RHSFirstPart = (velocity * velocity) / (gravity * deltaX);
        double RHSSecondPart = Math.Sqrt(
                (((velocity * velocity) * ((velocity * velocity) - (2 * gravity * deltaY)) 
                                / (gravity * gravity * deltaX * deltaX))
                                                - 1));
        double tanθ = RHSFirstPart - RHSSecondPart;
        double angle = Math.Atan2(tanθ, 1) * Mathf.Rad2Deg;
        if (angle < 0) return angle;
        return -angle;
    }
    
    Vector3 Mod(Vector3 Vec)
    {
        if (Vec.x < 0) Vec.x -= 2 * Vec.x;
        if (Vec.y < 0) Vec.y -= 2 * Vec.y;
        if (Vec.z < 0) Vec.z -= 2 * Vec.z;
        Vec.x /= 10;
        Vec.y /= 10;
        Vec.z /= 10;
        return Vec;
    }
}

Answer 1

Ok, as I can see, your implementation of formula from StackExchange is right, but you have to remember two things:

1. In unity there is a 3D world, so horizontal distance is not just pos1.x - pos2.x, but Mathf.Sqrt( deltaX * deltaX + deltaZ * deltaZ ), where deltaX = targetPositionMod.x - currentPosition.x and deltaZ = targetPositionMod.z - currentPosition.z
2. In computer implementation you have no 100% accuracy of math, so some problems can appear because of computational accuracy. And it can have affect on big distances. You can try to use double instead of float or find another implementation for arctangent function (I think, this can really help). But try this (second) advice only if first didn't help. It's harder to implement and it slows computations a bit.

Answer 2

Algorithm:

Step 1: Set up a function that calculates the appropriate solution of a quadratic equation

a*x^2 + b*x + c = 0

double quadratic_root(a,b,c){
   D = b^2 - 4*a*c
   return ( - b - Math.Sqrt(D) ) / (2 * a) 
}

Step 2: Input

current.x
current.y
current.z

target.x
target.y
target.z

velocity

gravity

Step 3: Calculate coefficients of the quadratic polynomial:

dist = Math.Sqrt( (target.x - current.x)^2 + (target.y - current.y)^2 )
a = gravity * dist^2 / (2 * velocity^2)
b = -dist
c = target.z - current.z + a

Step 4:

theta = Math.Atan2( quadratic_root(a,b,c), 1 )

Calculation behind the algorithm. You are in three space. The current position has coordinates

x = current.x
y = current.y
z = current.z

and the target has coordinates

x = target.x
y = target.y
z = target.z

Assume the angle between the initial velocity and the horizontal plane is theta. The magnitude of the projection of the distance between the current position and the target onto the horizontal $x,y-$plane is

dist = sqrt( (target.x - current.x)^2 - (target.y - current.y)^2 )

You are given the velocity magnitude velocity. Then, the speed with which the shadow (i.e. the orthogonal projection) of the arrow moves along the horizontal line between the source and the target is the magnitude of the shadow (i.e. the orthogonal projection) of the actual velocity

velocity * cos(theta)

The vertical speed of the arrow is then

velocity * sin(theta)

So the motion along dist follows the equation

dist = time * velocity * cos(theta)

and hence

time = dist / (velocity * cos(theta))

In the vertical direction, the motions is described by the equation

z = current.z + time * velocity * sin(theta) - time^2 * gravity / 2 

You are interested in the time for which the arrow hits the target, which has vertical coordinate target.z, so

target.z = current.z + time * velocity * sin(theta) - time^2 * gravity / 2 

The equation can be written as:

0 =  - (target.z - current.z) + time * velocity * sin(theta) - time^2 * gravity / 2 

We already know that

time = dist / (velocity * cos(theta))

so

0 = - (target.z - current.z) + dist * velocity * sin(theta) / (velocity * cos(theta)) - dist^2 * gravity / ( 2 * (velocity * cos(theta))^2 )  

which can be slightly simplified to

0 = - (target.z - current.z) + dist * sin(theta) / cos(theta) - gravity * dist^2 / ( 2 * (velocity * cos(theta))^2 )  

Because 1/( cos(theta)^2 ) = 1 + ( tan(theta) )^2 we obtain the quadratic in tan(theta) equation

a * ( tan(theta) )^2 + b * tan(theta) + c = 0

where

a = gravity * dist^2 / (2 * velocity^2)
b = - dist
c = target.z - current.z + a