How to use the poliastro.twobody.propagation.propagate function in poliastro

datetimeFormat = '%Y-%m-%d %H:%M:%S.%f'
    diff =datetime.strptime(str(f_time), datetimeFormat)\
        - datetime.strptime(str(time), datetimeFormat)    
    print("Time difference:")
    print(diff)
    tofs=TimeDelta(f_time-time)
    tofs = TimeDelta(np.linspace(0, 31* u.day, num=1))
    
 

    #print("tofs:")
    #print(tofs)
    #print("ie:")
    #print(initial.epoch.iso)
    rr =propagate(
        initial,
        tofs,
        method=cowell,
        ad=atmospheric_drag,
        R=R,
        C_D=C_D,
        A=A,
        m=m,
        H0=H0,
        rho0=rho0,
    )
  
    print("")
    print("Positions and velocity vectors are:")
    #print(str(rr.x))
    #print([float(s) for s in re.findall(r'-?\d+\.?\d*',str(rr.x))])

tofs,
            method=cowell,
            ad=a_d_1,
            #k=k,
            J2=J2,
            R=R,
            C_D=C_D,
            A=A,
            m=m,
            H0=H0,
            rho0=rho0
        )         

    elif(choice=='6'):

        rr =propagate(
            initial,
            tofs,
            method=cowell,
            ad=a_d_2,
            #k=k,
            J3=J3,
            R=R,
            C_D=C_D,
            A=A,
            m=m,
            H0=H0,
            rho0=rho0
        )          

    else:
        pass

Returns
        -------
        coordinate list
        """
        cart_cords = []  # type: List[float]

        h = (self.end_epoch - self.orbits[i][2]).to(u.second) / self.orbits[i][1]

        # Get rounding factor given the relative tolerance
        rf = 0
        while rtol < 1:
            rtol *= 10
            rf += 1

        for k in range(self.orbits[i][1] + 2):
            position = propagate(self.orbits[i][0], TimeDelta(k * h), rtol=rtol)

            cords = position.represent_as(CartesianRepresentation).xyz.to(u.meter).value
            cords = np.insert(cords, 0, h.value * k, axis=0)

            # flatten list
            cart_cords += list(map(lambda x: round(x[0], rf), cords.tolist()))

        return cart_cords

pass
    else:
        J3=input("Enter J3 constant:")
        R=input("Enter radius of earth in km:")       


    datetimeFormat = '%Y-%m-%d %H:%M:%S.%f'
    diff =datetime.strptime(str(f_time), datetimeFormat)\
        - datetime.strptime(str(time), datetimeFormat)    
    print("Time difference:")
    print(diff)
    tofs=TimeDelta(f_time-time)
    tofs = TimeDelta(np.linspace(0, 31* u.day, num=1))
    
    
    rr =propagate(
        initial,
        tofs,
        method=cowell,
        ad=J3_perturbation,
        J3=J3, 
        R=R
        )


    print("")
    print("Positions and velocity vectors are:")
    #print(str(rr.x))
    #print([float(s) for s in re.findall(r'-?\d+\.?\d*',str(rr.x))])
    r=[[float(s) for s in re.findall(r'-?\d+\.?\d*',str(rr.x))][0],[float(s) for s in re.findall(r'-?\d+\.?\d*',str(rr.y))][0],[float(s) for s in re.findall(r'-?\d+\.?\d*',str(rr.z))][0]]* u.km
    v=[[float(s) for s in re.findall(r'-?\d+\.?\d*',str(rr.v_x))][0],[float(s) for s in re.findall(r'-?\d+\.?\d*',str(rr.v_y))][0],[float(s) for s in re.findall(r'-?\d+\.?\d*',str(rr.v_z))][0]]* u.km / u.s
    f_orbit= Orbit.from_vectors(Earth, r, v)

def _init_orbit_packet_cords_(self, i):
        """

        Parameters
        ----------
        i: int
            Index of referenced orbit
        """
        h = (self.end_epoch - self.orbits[i][2]).to(u.second) / self.orbits[i][1]

        for k in range(self.orbits[i][1] + 2):
            position = propagate(self.orbits[i][0], TimeDelta(k * h))

            cords = position.represent_as(CartesianRepresentation).xyz.to(u.meter).value
            cords = np.insert(cords, 0, h.value * k, axis=0)

            self.czml[i]["position"]["cartesian"] += list(
                map(lambda x: x[0], cords.tolist())
            )

# Compute time of flight for correct epoch
            M = nu_to_M(self.nu, self.ecc)
            new_M = nu_to_M(value, self.ecc)
            time_of_flight = Angle(new_M - M).wrap_at(360 * u.deg) / self.n

            return self.from_classical(self.attractor, p / (1.0 - ecc ** 2) * u.km,
                                       ecc * u.one, inc * u.rad, raan * u.rad,
                                       argp * u.rad, value,
                                       epoch=self.epoch + time_of_flight, plane=self._plane)
        else:
            if isinstance(value, time.Time) and not isinstance(value, time.TimeDelta):
                time_of_flight = value - self.epoch
            else:
                time_of_flight = time.TimeDelta(value)

            return propagate(self, time_of_flight, method=method, rtol=rtol, **kwargs)

else:
        k_third=input("Enter moon's gravity:")
        x=input("Multiply Moon's gravity by X times so that effect is visible,input X:")    

    datetimeFormat = '%Y-%m-%d %H:%M:%S.%f'
    diff =datetime.strptime(str(f_time), datetimeFormat)\
        - datetime.strptime(str(time), datetimeFormat)    
    print("Time difference:")
    print(diff)
    tofs=TimeDelta(f_time-time)
    tofs = TimeDelta(np.linspace(0, 31* u.day, num=1))
    
    

    # multiply Moon gravity by x so that effect is visible :)
    rr = propagate(
        initial,
        tofs,
        method=cowell,
        rtol=1e-6,
        ad=third_body,
        k_third=x *k_third,
        third_body=body_r,
    )

    print("")
    print("Positions and velocity vectors are:")
    #print(str(rr.x))
    #print([float(s) for s in re.findall(r'-?\d+\.?\d*',str(rr.x))])
    r=[[float(s) for s in re.findall(r'-?\d+\.?\d*',str(rr.x))][0],[float(s) for s in re.findall(r'-?\d+\.?\d*',str(rr.y))][0],[float(s) for s in re.findall(r'-?\d+\.?\d*',str(rr.z))][0]]* u.km
    v=[[float(s) for s in re.findall(r'-?\d+\.?\d*',str(rr.v_x))][0],[float(s) for s in re.findall(r'-?\d+\.?\d*',str(rr.v_y))][0],[float(s) for s in re.findall(r'-?\d+\.?\d*',str(rr.v_z))][0]]* u.km / u.s
    f_orbit= Orbit.from_vectors(Earth, r, v)

if(check=='1'):
        pass
    else:
        J2=input("Enter J2 constant")
        R=input("Enter radius of earth in km")        



    datetimeFormat = '%Y-%m-%d %H:%M:%S.%f'
    diff =datetime.strptime(str(f_time), datetimeFormat)\
        - datetime.strptime(str(time), datetimeFormat)    
    print("Time difference:")
    print(diff)
    tofs=TimeDelta(f_time-time)
    tofs = TimeDelta(np.linspace(0, 31* u.day, num=1))
    rr =propagate(
        initial,
        tofs,
        method=cowell,
        ad=J2_perturbation,
        J2=J2, 
        R=R
        )


    print("")
    print("Positions and velocity vectors are:")
    #print(str(rr.x))
    #print([float(s) for s in re.findall(r'-?\d+\.?\d*',str(rr.x))])
    r=[[float(s) for s in re.findall(r'-?\d+\.?\d*',str(rr.x))][0],[float(s) for s in re.findall(r'-?\d+\.?\d*',str(rr.y))][0],[float(s) for s in re.findall(r'-?\d+\.?\d*',str(rr.z))][0]]* u.km
    v=[[float(s) for s in re.findall(r'-?\d+\.?\d*',str(rr.v_x))][0],[float(s) for s in re.findall(r'-?\d+\.?\d*',str(rr.v_y))][0],[float(s) for s in re.findall(r'-?\d+\.?\d*',str(rr.v_z))][0]]* u.km / u.s
    f_orbit= Orbit.from_vectors(Earth, r, v)