Using Latitude, Longitude, Bearing and Distance to get resulting coordinates giving inconsistent elliptical behaviour

I have the following inputs

1. Latitude and Longitude in decimal notation
2. Bearing (in degrees), distance in km

I want to find out the latitude and longitude of the point which lies at the specified distance and bearing from the original point (I'm looking using distances in the range of 0-1000 meters, and the location is pretty close to the equator).

I used the following function

def calculate_destination_point(lat, lon, distance, bearing):
    """Input: latitude, longitude, distance in kilometers, bearing of line
    Output: lat, long of the points at specified distance and bearing.
    """
    # Constants
    radius = 6371  # Radius of the Earth in kilometers

    # Convert distance to radians
    distance /= radius

    # Convert bearing to radians
    bearing = math.radians(bearing)

    # Convert latitude and longitude to radians
    lat_rad = math.radians(lat)
    lon_rad = math.radians(lon)

    # Calculate destination point latitude
    dest_lat = math.asin(math.sin(lat_rad) * math.cos(distance) +
                         math.cos(lat_rad) * math.sin(distance) * math.cos(bearing))

    # Calculate destination point longitude
    dest_lon = lon_rad + math.atan2(math.sin(bearing) * math.sin(distance) * math.cos(lat_rad),
                                    math.cos(distance) - math.sin(lat_rad) * math.sin(dest_lat))

    # Convert latitude and longitude back to degrees
    dest_lat = math.degrees(dest_lat)
    dest_lon = math.degrees(dest_lon)

    return dest_lat, dest_lon

For inputs

bearing: 45
distance(km): 0.02
Source: 78.720249,9.944548

My output is

Destination Lat Lon: 78.72037618257347,9.945198229935736

I'm able to validate the distance with haversine's formula and the bearing using

def calculate_bearing(lat1, lon1, lat2, lon2):
    lat1_rad = math.radians(lat1)
    lat2_rad = math.radians(lat2)
    lon_diff_rad = math.radians(lon2 - lon1)

    y = math.sin(lon_diff_rad) * math.cos(lat2_rad)
    x = math.cos(lat1_rad) * math.sin(lat2_rad) - math.sin(lat1_rad) * math.cos(lat2_rad) * math.cos(lon_diff_rad)

    bearing_rad = math.atan2(y, x)
    bearing_deg = math.degrees(bearing_rad)
    bearing_deg = (bearing_deg + 360) % 360

    return bearing_deg

But when I convert these coordinates to a KML file and plot it, the following happened For bearings 0 and 180, the measured distance in google earth showed up as the correct 0.02 value, but for bearings 90 and 180, it showed up 0.10, and for 45, it was 0.07.

There seems to be some sort of elliptical behaviour around a single point. As the bearing increases from 0 to 180, the coordinate returned is further away from the source point.

Here's a snippet of the code I'm using to get the elliptical behaviour. The main thing to note is that only the bearing is changing.

    lat1, lon1 = 78.720249,9.944548
    distance = 0.02
    kml = simplekml.Kml()
    for bearing in range(0,360,15):
        coords = calculate_destination_point(lat1,lon1,distance,bearing)
        kml.newlinestring(name=f"Bearing:{bearing}",coords=[(lat1,lon1),coords])
    
    kml.save('../nalla/kml_practice/kml_dummies/elliptical.kml')

I've attached a screenshot of the behaviour.

If this is normal behaviour then how do I correct it? And if it isn't could someone find the fault?

What did I try: StackOverflow showed me a post about using GeographicLib, but that me an even bigger ellipse.

Part of me is wondering if I have the correct results but incorrect processing of it.

Here is the output from the GeographicLib-based functions. As you would observe, the calculated distance is the same, but when plotting it using KML, I'm getting an ellipse. Ellipses Endpoint Lat: 78.72058074864637 Endpoint Lon: 9.944548 Endpoint Bearing: 0 Endpoint Distance: 0.03704000000088827

Endpoint Lat: 78.72056944426299 Endpoint Lon: 9.94498687182839 Endpoint Bearing: 15 Endpoint Distance: 0.03704000000059051

Endpoint Lat: 78.7205363015523 Endpoint Lon: 9.945395832808659 Endpoint Bearing: 30 Endpoint Distance: 0.037040000000551486

Endpoint Lat: 78.72048357931116 Endpoint Lon: 9.94574701112629 Endpoint Bearing: 45 Endpoint Distance: 0.037040000000121844

---

Finally: Turns out the error was purely because of the ordering of lat longs. After looking at the code with a fresh mind, I was able to correct it with the help of the answer marked correct, and the comment from @Ture Palsson! Appreciate everyone's earnest reply and help!

The issue is here:

kml.newlinestring(name=f"Bearing:{bearing}",coords=[(lat1,lon1),coords])

The SimpleKML docs say that the order of the coordinates is lon, lat:

The coordinates of the feature, accepts list of tuples.

A tuple represents a coordinate in the order longitude then latitude. The tuple has the option of specifying a height. If no height is given, it defaults to zero. A point feature has just one point, therefore a list with one tuple is given.

Examples:

• No height: [(1.0, 1.0), (2.0, 1.0)]
• Height: [(1.0, 1.0, 50.0), (2.0, 1.0, 10.0)]
• Point: [(1.0, 1.0)] # longitude, latitude
• Point + height: [(1.0, 1.0, 100)] # longitude, latitude, height of 100m