The Collapse of Early Neolithic Settlements in the Southern Levant (Israel)

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

1. --- DATA (sama kuin aiemmin) ---

years = np.array([7250,7150,7050,6950,6850,6750,6650,6550,6450,6350,6250,6150,6050,5950,5850,5750]) houses = np.array([126,157,188,219,250,283,346,409,472,535,600,637,674,711,748,787])

ref_years = np.array([7250,6750,6250,5750]) families_avg = np.array([97.5,215.5,464,612.5]) ref_hectares = np.array([195,441,928,1225]) average_family_size = 6

1. Population

houses_at_ref = houses0,5,10,15 ratio_fam = families_avg / houses_at_ref interp_ratio_fam = np.interp(years, ref_years, ratio_fam) population = houses * interp_ratio_fam * average_family_size

1. Radius (km)

ratio_ha = ref_hectares / houses_at_ref interp_ratio_ha = np.interp(years, ref_years, ratio_ha) total_hectares = houses * interp_ratio_ha radius_km = np.sqrt(total_hectares * 10000 / np.pi) / 1000 # metri → km

time_elapsed = years[0] - years

1. Määritellään logistinen kasvufunktio
2. P(t) = K / (1 + ((K - P0) / P0) * e^(-rt))

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

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

1. Sovitetaan käyrä dataan.
2. Annetaan alkuarvaukset (p0) auttamaan laskentaa:
3. K=5000 (arvaus), r=0.002 (hidas kasvu), P0=population[0]

try:

   popt, pcov = curve_fit(logistic_growth, time_elapsed, population, 
                          p0=[5000, 0.002, population[0]], 
                          bounds=([max(population), 0, 0], [20000, 0.05, 2000]))
   
   K_estimated = popt[0] # Ensimmäinen parametri on K
   r_overall = popt[1]
   
   k_text = f"Estimated K $\\approx$ {int(K_estimated)}"
   model_success = True

except:

   K_estimated = max(population)
   k_text = "Cannot estimete K, rapid growth"
   model_success = False

1. print(K_estimated)

1. quit(-1)

1. Logistic fit

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

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

t_elapsed = years[0] - years

1. popt, pcov = curve_fit(logistic, t_elapsed, population,
2. p0=[12500, 0.0027, population[0]],
3. bounds=([7000, 0.0005, 300], [35000, 0.01, 3000]),
4. maxfev=10000)

K, r0, P0 = popt

model = logistic(t_elapsed, *popt)

1. ==================== KUVAAJA ====================

fig, ax = plt.subplots(figsize=(16, 10), dpi=120)

1. 1. Millimetripaperimainen tausta

ax.set_facecolor('#f8f8ff') ax.grid(True, which='minor', color='#c0d6e4', linewidth=0.5, alpha=0.7) ax.grid(True, which='major', color='#88aacc', linewidth=0.9, alpha=0.9)

ax.set_axisbelow(True)

1. Asetetaan tiheä minor-ruudukko (50 vuoden ja 500 hengen välein)

ax.set_xticks(np.arange(7300, 5600, -50), minor=True) ax.set_yticks(np.arange(0, 40000, 500), minor=True) ax.set_xticks(np.arange(7300, 5600, -200)) ax.set_yticks(np.arange(0, 40000, 2000))

1. 2. Pääkäyrät

ax.plot(years, population, 'o', color='#1a9850', markersize=10, markeredgecolor='darkgreen', markeredgewidth=1.5,

       label='Estimated population', zorder=10)

ax.plot(years, model, '--', color='#006400', linewidth=4.5, label='Logistic model', zorder=9) ax.axhline(K_estimated, color='#cc0000', linestyle=':', linewidth=4, alpha=0.9,

          label=f'Carrying capacity K ≈ {int(K_estimated):,}')

1. 3. Viljelyalueen säde (toissijainen y-akseli)

ax2 = ax.twinx() ax2.fill_between(years, 0, radius_km, color='#8c564b', alpha=0.12, zorder=1) ax2.plot(years, radius_km, color='#8c564b', linewidth=4, label='Cultivated radius (km)') ax2.set_ylabel('Cultivated radius (km)', color='#8c564b', fontsize=16, fontweight='bold') ax2.tick_params(axis='y', labelcolor='#8c564b', labelsize=13) ax2.set_ylim(0, 4.2)

1. 4. Paikalliset annotaatiot (kuten aiemmin, mutta nätimmässä laatikossa)

targets = [7000, 6500, 6000] for yr in targets:

   t = years[0] - yr
   P = logistic(t, *popt)
   r_local = r0 * (1 - P/K_estimated)
   dbl = np.log(2)/r_local if r_local>1e-6 else 99999

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

1. 5. Romahtamisvyöhyke ja teksti

collapse_left = 5900 collapse_right = 5700 ax.axvspan(collapse_right, collapse_left, alpha=0.25, color='#ff9999', zorder=0) ax.text(5810, K_estimated*0.78*0.5, "Settlement\nabandonment\n~5850–5750 BC", ha='center', va='center',

       fontsize=18, fontweight='bold', color='#990000',
       bbox=dict(boxstyle="round,pad=0.8", facecolor="#ffcccc", edgecolor="#990000", lw=2))

1. 6. Otsikko ja legendat

ax.set_title("’Ayn Ghazal (Jordan) – Neolithic demographic boom and collapse\n"

            "Logistic growth + cultivated radius + local doubling time", 
            fontsize=22, pad=25, fontweight='bold')

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

1. Yhdistetty legenda

lines1, labels1 = ax.get_legend_handles_labels() lines2, labels2 = ax2.get_legend_handles_labels() ax.legend(lines1 + lines2, labels1 + labels2, loc='upper left', fontsize=14, framealpha=1, fancybox=True)

plt.tight_layout() plt.savefig("ayn_ghazal_population_change_1_1_1_1.png")

plt.show()

1. Konsoliin tarkat luvut

print(f"K = {K:,.0f} | r₀ = {r0*100:.4f}%/yr") for yr in targets:

   t = years[0]-yr
   P = logistic(t,*popt)
   print(f"{yr} BC → {int(P):,} people | r(t) = {r0*(1-P/K_estimated)*100:.3f}%/yr | doubling {np.log(2)/(r0*(1-P/K_estimated)):.0f} yr")