def calculate_pv_system_production(
    pv_system: SimplifiedPVSystem, solar_inputs: SolarInputs
) -> np.ndarray:
    """
    Returns AC power output in kW using the PVWatts model.
    """
    # 1. Location
    location = pvlib.location.Location(
        latitude=pv_system.coordinates.latitude,
        longitude=pv_system.coordinates.longitude,
        tz=str(pv_system.timezone),
        altitude=pv_system.coordinates.altitude,
    )

    # 2. System
    system = pvlib.pvsystem.PVSystem(  # todo review all of these parameters with dict
        surface_tilt=pv_system.surface_tilt,
        surface_azimuth=pv_system.surface_azimuth,
        module_parameters={
            "pdc0": pv_system.solar_power_kwp * 1000,  # W peak DC
            "gamma_pdc": pv_system.temp_coefficient,
        },
        inverter_parameters={
            "pdc0": pv_system.solar_power_kwp * 1000,
            "eta_inv_nom": pv_system.inverter_efficiency,
        },
        temperature_model_parameters=pvlib.temperature.TEMPERATURE_MODEL_PARAMETERS[
            "sapm"
        ]["open_rack_glass_glass"],
    )

    # 3. ModelChain
    mc = pvlib.modelchain.ModelChain(
        system,
        location,
        dc_model="pvwatts",
        ac_model="pvwatts",
        aoi_model="physical",
        spectral_model="no_loss",
        losses_model="pvwatts",
    )

    # 4. Run model — pvlib computes sun position internally from location + times
    weather_data = {
        "dni": solar_inputs.dni_wm2,
        "dhi": solar_inputs.dhi_wm2,
        "ghi": solar_inputs.ghi_wm2,
    }
    if solar_inputs.outdoor_air_temp_c is not None:
        weather_data["temp_air"] = solar_inputs.outdoor_air_temp_c

    solar_inputs_pvlib_format = pd.DataFrame(
        index=solar_inputs.ts,
        data=weather_data,
    )

    mc.run_model(solar_inputs_pvlib_format)

    # 5. Return AC power in kW
    return (mc.results.ac / 1000).to_numpy()