Adjusted for stare mode of gianob

        # ======================== BARYCENTRIC VELOCITY EXTRACTION ========================

        # Initialize BERV (Barycentric Earth Radial Velocity) values

        self.berv = -1      # BERV at observation time

        self.bervmax = -1   # Maximum BERV during exposure

        self.berv = None      # BERV at observation time

        self.bervmax = None   # Maximum BERV during exposure

        # Extract BERV from header if available

        if yaml_file["BERV"] != "None":

            self.berv = self.header[yaml_file["BERV"]]

import sys

import os

import numpy as np

from PyAstronomy import pyasl

import pickle

import matplotlib

from astropy.time import Time

import scipy.io as sio

from pathlib import Path

import numpy as np

from astropy.time import Time

from astropy.coordinates import SkyCoord, EarthLocation, solar_system_ephemeris, get_body_barycentric_posvel

import astropy.units as u

def bjd_calculator(mjd, ra, dec, lon, lat, alt):

    """

    Calculate BJD and barycentric correction

    Parameters:

    -----------

    mjd : float or array - Modified Julian Date

    ra : float or array - Right Ascension (degrees)

    dec : float or array - Declination (degrees)

    lon : float - Longitude (degrees, East positive)

    lat : float - Latitude (degrees)

    alt : float - Altitude (meters)

    Returns:

    --------

    bjd : float or array - Barycentric Julian Date

    bcorr : float or array - Barycentric correction (days)

    """

    rjd = mjd + 0.5 # reduced julian date

    jd = rjd + 2400000 # UTC julian date

    time = Time(jd, format='jd', scale='utc')

    location = EarthLocation.from_geodetic(lon * u.deg, lat * u.deg, alt * u.m)

    target = SkyCoord(ra * u.deg, dec * u.deg, frame='icrs')

    ltt_bary = time.light_travel_time(target, 'barycentric', location=location)

    bjd = time.tdb + ltt_bary

    return bjd.jd, ltt_bary.to(u.day).valueù

def helcorr_velocity(mjd, ra, dec, lon, lat, alt):

    """

    Calculate barycentric velocity correction like IDL's helcorr

    Parameters:

    -----------

    mjd : float or array - Modified Julian Date

    ra : float - Right Ascension (degrees)

    dec : float - Declination (degrees)

    lon : float - Longitude (degrees, East positive)

    lat : float - Latitude (degrees)

    alt : float - Altitude (meters)

    Returns:

    --------

    corr : float or array - Barycentric velocity correction in km/s

        This is the velocity to ADD to observed radial velocity

    hjd : float or array - Heliocentric Julian Date (full JD)

    vbar : float or array - Barycentric velocity in km/s

    """

    rjd = mjd + 0.5  # reduced julian date

    jd = rjd + 2400000  # UTC julian date

    # Create Time object

    time = Time(jd, format='jd', scale='utc')

    # Observatory location

    location = EarthLocation.from_geodetic(lon * u.deg, lat * u.deg, alt * u.m)

    # Target coordinates

    target = SkyCoord(ra * u.deg, dec * u.deg, frame='icrs')

    # Calculate HJD (ltt = light time travel)

    ltt_helio = time.light_travel_time(target, 'heliocentric', location=location)

    hjd = (time.utc + ltt_helio).jd

    # Get Earth's barycentric velocity

    with solar_system_ephemeris.set('builtin'):

        earth_pv = get_body_barycentric_posvel('earth', time)

    # Velocity in km/s

    v_earth = earth_pv[1].xyz.to(u.km / u.s).value

    # Project velocity onto line of sight

    ra_rad = np.radians(ra)

    dec_rad = np.radians(dec)

    vbar = (v_earth[0] * np.cos(dec_rad) * np.cos(ra_rad) +

            v_earth[1] * np.cos(dec_rad) * np.sin(ra_rad) +

            v_earth[2] * np.sin(dec_rad))

    # Calculate diurnal velocity (Earth's rotation)

    # This is the contribution from observatory location

    # Observer's velocity in Earth frame

    v_obs = location.get_gcrs_posvel(time)[1]

    # Project onto target direction

    vdiurnal = (v_obs.x.to(u.km / u.s).value * np.cos(dec_rad) * np.cos(ra_rad) +

                v_obs.y.to(u.km / u.s).value * np.cos(dec_rad) * np.sin(ra_rad) +

                v_obs.z.to(u.km / u.s).value * np.sin(dec_rad))

    # Total correction (this is what helcorr returns as 'corr')

    corr = vdiurnal + vbar

    return corr, hjd, vbar

def main():

    path_default = sys.argv[9]

    sys.path.append(path_default)

    os.chdir(path_default)

            corri, hjdi = pyasl.helcorr(

                long, lat, alt, ra[i], dec[i], rjd[i] + 2.4e6, debug=False

            )

        if berv[i] > 0:

            print(f"Order {i}: berv calc = {corri}, berv header = {berv[i]}, bervmax = {bervmax[i]}")

            corri2, hjdi2, vbar2 = helcorr_velocity(mjd[i], ra[i], dec[i], long, lat, alt)

            # print(f"Compare old: {corri} vs new {corri2} for mjd {mjd[i]}. Difference new - old {corri2 - corri}")

            if berv[i] is not None:

                print(f"Order {i}: berv calc1 = {corri}, berv header = {berv[i]}, bervmax = {bervmax[i]}")

                corri = berv[i]

            if corr_zero:

                corri = 0

        with open(str(Path(path, pref + "phase.pkl")), "wb") as fof:

            pickle.dump(phase_dictionary, fof)

        print("From header 4/4 OK")

if __name__ == "__main__":

    main()