import numpy as np import matplotlib.pyplot as plt from scipy.optimize import curve_fit

1. --- FULL PHASE DATA (same as before) ---

years_bc = np.array([7100,7050,7000,6950,6900,6850,6800,

                    6750,6700,6650,6600,6550,6500,
                    6450,6400,6350,6300,6250,
                    6200,6150,6100,6050,6000,5950])

population = np.array([150,250,350,450,550,650,700,

                      720,750,740,730,720,700,
                      650,550,450,350,300,
                      250,200,150,120,100,50])

1. --- ONLY GROWTH PHASE FOR LOGISTIC FIT (7100–6500 BC) ---

mask_growth = years_bc >= 6500 years_growth = years_bc[mask_growth] pop_growth = population[mask_growth]

time_elapsed_growth = years_growth[0] - years_growth # t=0 at 7100 BC

def logistic(t, K, r, P0):

   return K / (1 + ((K - P0)/P0) * np.exp(-r * t))

1. Fit only to growth phase → realistic K

popt, pcov = curve_fit(logistic, time_elapsed_growth, pop_growth,

                      p0=[1100, 0.0018, 150],
                      bounds=([800, 0.0005, 50], [2000, 0.01, 400]),
                      maxfev=5000)

K_est, r_intr, P0_est = popt

1. Extend model to full timeline for visualization

time_full = years_bc[0] - years_bc model_full = logistic(time_full, *popt)

1. --- PLOT ---

fig, ax = plt.subplots(figsize=(16, 10), dpi=130) ax.set_facecolor('#f8f8ff') ax.grid(True, which='minor', color='#c0d6e4', lw=0.5, alpha=0.7) ax.grid(True, which='major', color='#88aacc', lw=0.9, alpha=0.9) ax.set_axisbelow(True)

ax.set_xticks(np.arange(7200, 5900, -50), minor=True) ax.set_yticks(np.arange(0, 1600, 100), minor=True) ax.set_xticks(np.arange(7200, 5900, -200)) ax.set_yticks(np.arange(0, 1600, 200))

1. All data points

ax.plot(years_bc, population, 'o', color='#1a9850', ms=11, mew=2, mec='darkgreen',

       label='Estimated Population (Phase Data)', zorder=10)

1. Logistic model (dashed, only growth phase solid for emphasis)

ax.plot(years_growth, logistic(time_elapsed_growth, *popt),

       '-', color='#006400', lw=5, label='Logistic Model (fitted 7100–6500 BC)', zorder=9)

ax.plot(years_bc[years_bc < 6500], model_full[years_bc < 6500],

       '--', color='#006400', lw=4, alpha=0.7, zorder=9)

ax.axhline(K_est, color='#cc0000', linestyle=':', linewidth=4, alpha=0.9,

          label=f'Carrying Capacity K ≈ {int(K_est):,} (growth-phase estimate)')

1. Annotations only in growth + peak phase

targets = [7000, 6700, 6500] # Now 6500 BC is near the real peak! for yr in targets:

   t = years_bc[0] - yr
   P_model = logistic(t, *popt)
   P_data  = population[np.abs(years_bc - yr).argmin()]
   r_local = r_intr * (1 - P_model / K_est)
   dbl = np.log(2) / r_local if r_local > 1e-6 else 99999

   txt = f"{yr} BC\nPop ≈ {int(round(P_data)):,}\nr(t) ≈ {r_local*100:.3f} %/yr\nDoubling ≈ {dbl:.0f} yr"
   ax.annotate(txt, xy=(yr, P_data), xytext=(35, 45), textcoords='offset points',
               bbox=dict(boxstyle="round,pad=0.7", fc="#e8f5e8", ec="#004d00", lw=2),
               fontsize=13.5, fontweight='bold', color='#004d00',
               arrowprops=dict(arrowstyle='->', color='#004d00', lw=2.5))

1. Decline zone

ax.axvspan(6300, 5950, alpha=0.25, color='#ff9999', zorder=0) ax.text(6120, K_est*0.75, "Population decline\n& abandonment\n~6300–5950 BC",

       ha='center', va='center', fontsize=17, fontweight='bold', color='#990000',
       bbox=dict(boxstyle="round,pad=0.9", fc="#ffcccc", ec="#990000", lw=2.5))

1. Title & labels

ax.set_title("Çatalhöyük – estimated Neolithic Population Dynamics\n"

            "Logistic model fitted only to growth phase (7100–6500 BC) – Kuijt & Marciniak 2024",
            fontsize=20, pad=25, fontweight='bold')

ax.set_xlabel("Year BC", fontsize=20, fontweight='bold') ax.set_ylabel("Population", fontsize=20, fontweight='bold', color='#1a9850') ax.tick_params(axis='y', labelcolor='#1a9850', labelsize=18) ax.tick_params(axis='x', labelsize=18) ax.invert_xaxis() ax.set_ylim(0, 1350)

1. Info box: why we fit only growth phase

ax.text(0.02, 0.98, "Note: Logistic model fitted only to Early+Middle phase\n"

                   "Late/Final decline driven by environmental factors,\n"
                   "not simple density-dependent growth",
       transform=ax.transAxes, fontsize=14, verticalalignment='top',
       bbox=dict(boxstyle="round", fc="wheat", ec="orange", alpha=0.9))

ax.legend(loc='upper left', fontsize=16, framealpha=1, fancybox=True)

plt.tight_layout()

plt.savefig("catal_hoyuk_estimated_population_1_1_1_1.png") plt.show()

1. Console output

print("Çatalhöyük – Logistic fit ONLY to growth phase (7100–6500 BC)") print(f"K ≈ {K_est:,.0f} | r₀ ≈ {r_intr*100:.4f}%/yr") for yr in targets:

   t = years_bc[0] - yr
   P = logistic(t, *popt)
   r_loc = r_intr * (1 - P / K_est)
   dt = np.log(2) / r_loc
   print(f"{yr} BC → Pop ≈ {int(round(P)):,} | r(t) ≈ {r_loc*100:.3f}%/yr | doubling ≈ {dt:.0f} yr")