---
jupytext:
  formats: md:myst
  text_representation:
    extension: .md
    format_name: myst
kernelspec:
  display_name: Python 3
  language: python
  name: python3
---

(har-joaca-te)=
# {fa}`rocket` Calculează tu: stereograma, înclinarea aparentă, cele trei puncte

Trei calcule scurte din acest capitol, pe care le poți rula și modifica. **Schimbă numerele din prima celulă a fiecărui exemplu** și vezi ce se întâmplă.

## Un plan pe rețeaua stereografică

Planul este dat prin *sens de înclinare/înclinare*. Codul desenează urma lui ciclografică și polul, în proiecție stereografică pe emisfera inferioară, cu relația {eq}`eq-har-stereo-raza`.

```{code-cell} ipython3
sens, inclinare = 120, 35     # poziția planului: sens de înclinare / înclinare [grade]
```

```{code-cell} ipython3
:tags: [hide-input]
import numpy as np
import matplotlib.pyplot as plt

def proiectie(azimut, plonjare):
    """linie (azimut, plonjare) -> (x, y) pe rețeaua de rază 1; nordul în sus"""
    r = np.tan(np.radians(90 - plonjare) / 2)
    a = np.radians(azimut)
    return r * np.sin(a), r * np.cos(a)

def urma_ciclografica(sens, inclinare, n=181):
    """liniile conținute în plan: unghiul t se măsoară în plan, de la direcție"""
    t = np.radians(np.linspace(0, 180, n))
    d, s = np.radians(inclinare), np.radians(sens)
    # vector în plan: combinație între direcție (orizontală) și linia de cea mai mare pantă
    e = np.cos(t) * np.sin(s - np.pi / 2) + np.sin(t) * np.cos(d) * np.sin(s)   # est
    n_ = np.cos(t) * np.cos(s - np.pi / 2) + np.sin(t) * np.cos(d) * np.cos(s)  # nord
    jos = np.sin(t) * np.sin(d)                                                 # în jos
    az = np.degrees(np.arctan2(e, n_)) % 360
    pl = np.degrees(np.arcsin(np.clip(jos, 0, 1)))
    return proiectie(az, pl)

fig, ax = plt.subplots(figsize=(5, 5))
u = np.linspace(0, 2 * np.pi, 361)
ax.plot(np.sin(u), np.cos(u), "k", lw=1)
for az in range(0, 360, 30):
    ax.plot(*proiectie(np.array([az, az]), np.array([0, 3])), "k", lw=0.8)
for txt, az in (("N", 0), ("E", 90), ("S", 180), ("V", 270)):
    ax.text(*(1.09 * np.array([np.sin(np.radians(az)), np.cos(np.radians(az))])), txt, ha="center", va="center")

ax.plot(*urma_ciclografica(sens, inclinare), color="tab:red", lw=2, label=f"planul {sens:03d}/{inclinare}")
pol_az, pol_pl = (sens + 180) % 360, 90 - inclinare
ax.plot(*proiectie(pol_az, pol_pl), "o", color="tab:red", label=f"polul: azimut {pol_az:03d}°, plonjare {pol_pl}°")
ax.plot(*proiectie(sens, inclinare), "s", color="tab:blue", label="linia de cea mai mare pantă")
ax.plot(0, 0, "k+")
ax.set_aspect("equal"); ax.axis("off")
ax.legend(loc="upper center", bbox_to_anchor=(0.5, -0.02), frameon=False)
plt.show()
print(f"polul se află la r = {np.tan(np.radians(inclinare) / 2):.3f} R de centru")
```

Încearcă `inclinare = 0` (plan orizontal: urma este cercul primitiv, polul este în centru) și `inclinare = 90` (plan vertical: urma este un diametru, polul este pe margine).

## Înclinarea aparentă

Cât de înclinat *pare* un strat într-o secțiune care nu este trasată pe sensul lui de înclinare? Relația {eq}`eq-har-inclinare-aparenta`.

```{code-cell} ipython3
alpha = 40        # înclinarea reală [grade]
```

```{code-cell} ipython3
:tags: [hide-input]
beta = np.arange(0, 91, 15)
ap = np.degrees(np.arctan(np.tan(np.radians(alpha)) * np.cos(np.radians(beta))))
print("unghiul dintre secțiune și sensul de înclinare  ->  înclinarea aparentă")
for b, a in zip(beta, ap):
    print(f"   {b:3d}°   ->   {a:5.1f}°")

b = np.linspace(0, 90, 91)
plt.figure(figsize=(6, 3.2))
plt.plot(b, np.degrees(np.arctan(np.tan(np.radians(alpha)) * np.cos(np.radians(b)))))
plt.xlabel("β – unghiul dintre traseul secțiunii și sensul de înclinare [°]")
plt.ylabel("înclinarea aparentă [°]"); plt.grid(alpha=0.3)
plt.show()
```

Până la β ≈ 30° eroarea este mică; dincolo de 60° stratul pare aproape de două ori mai puțin înclinat decât este.

## Problema celor trei puncte

Trei puncte în care a fost întâlnită aceeași suprafață (de exemplu baza unui strat, în trei foraje) definesc un plan. Codul îl găsește pentru orice trei puncte necoliniare; aplicația A2 din {ref}`sec-har-aplicatii` este același calcul, făcut grafic.

```{code-cell} ipython3
# (x spre est [m], y spre nord [m], cota suprafeței [m])
A = (100.0, 400.0, 400.0)
B = (500.0, 450.0, 250.0)
C = (350.0, 100.0, 250.0)
```

```{code-cell} ipython3
:tags: [hide-input]
p1, p2, p3 = map(np.array, (A, B, C))
n = np.cross(p2 - p1, p3 - p1)          # normala la plan
if abs(n[2]) < 1e-9:
    raise ValueError("puncte coliniare (în plan) sau plan vertical: schimbă datele")
if n[2] < 0:
    n = -n                              # normala orientată în sus
# planul coboară în sensul componentei orizontale a normalei
sens = np.degrees(np.arctan2(n[0], n[1])) % 360
inclinare = np.degrees(np.arctan(np.hypot(n[0], n[1]) / n[2]))
print(f"poziția planului: {sens:03.0f}/{inclinare:.0f}   (sens de înclinare / înclinare)")
print(f"direcția: {(sens - 90) % 180:03.0f}° – {(sens - 90) % 180 + 180:03.0f}°")
grad = np.tan(np.radians(inclinare))
print(f"izobate din 50 în 50 m: distanța orizontală dintre ele = {50 / grad:.0f} m")
```

Mută punctul C pe dreapta AB și rulează din nou: calculul refuză, pentru că trei puncte coliniare nu definesc un plan (vezi și exercițiul 8 al capitolului).