#!/usr/bin/env python3
"""Gera deta_50_1mb — distribuicao vertical LM=50 para topo PT=1 hPa.

Desenho (2026-08-06, ver memoria eta-tv/1hPa):
  - Zona S (L1-11, 1->87 hPa): Delta(ln p) afinando 0.70->0.27 (dz ~5.3->1.7 km),
    monotonico em dz => sem degrau de refletividade p/ ondas de gravidade; topo
    relativamente espesso (buffer da tampa rigida, complementa CODAMP/DDAMP).
  - Zona M (L12-37, 87->888 hPa): parametrica, sobe suave (razao <=1.09) ate
    ~42 hPa em 500 hPa e decai ate casar com a PBL; a subida lenta em 100-300 hPa
    preserva refinamento relativo na tropopausa/jato.
  - Zona B (L38-50, 888->1000 hPa): as 13 ultimas camadas do deta_50_25mb
    VERBATIM em dp (fisica de superficie/SFCDIF identica ao operacional).
Formato do arquivo: 1 registro unformatted little-endian = DETA(50) R4 +
LDUM(51) I4 (dummy, copiado do baseline; interp.f le e ignora).
Custo da tampa a 1 hPa com LM=50: media troposfera ~42 hPa vs ~28 no baseline
(400-600 hPa). Alternativa estrutural sem custo troposferico: LM=60 (PARMETA +
tabelas CO2 novas).
"""
import numpy as np, struct, sys, os

SUP = '/home/jorge/Eta_support_data/static/fix/deta_files'
base = os.path.join(SUP, 'deta_50_25mb')
out  = os.path.join(SUP, 'deta_50_1mb')

raw = open(base, 'rb').read()
nrec = struct.unpack('<i', raw[:4])[0]
d25  = np.frombuffer(raw[4:204], dtype='<f4').astype(float)
tail = raw[204:4+nrec]                      # LDUM(51) dummy
PT25, PS = 25., 1000.
dp25 = d25*(PS-PT25)

PT = 1.; COL = PS-PT
# Zona S
dlnp = np.array([0.70,0.58,0.49,0.44,0.40,0.37,0.34,0.31,0.29,0.28,0.27])
pS  = PT*np.exp(np.concatenate([[0], np.cumsum(dlnp)]))
dpS = np.diff(pS)
# Zona B
nbot = 13
dpB  = dp25[-nbot:]
# Zona M
nmid = 50-len(dpS)-nbot
span = (PS-dpB.sum())-pS[-1]
basel = np.linspace(dpS[-1]*1.05, dpB[0]/0.94, nmid)
x = np.linspace(0, 1, nmid)
shape = np.sin(np.pi*x)**1.5
dpM = basel + (span-basel.sum())/shape.sum()*shape

dp   = np.concatenate([dpS, dpM, dpB])
deta = (dp/COL).astype(np.float32)          # residuo ~1e-7: interp.f forca eta(LM+1)=1.

rec = deta.tobytes() + tail
open(out, 'wb').write(struct.pack('<i', len(rec)) + rec + struct.pack('<i', len(rec)))

pint = np.concatenate([[PT], PT+np.cumsum(dp)])
print(f"gravado {out}: LM=50 PT=1hPa  sum={deta.sum(dtype=np.float32):.7f}")
print(f"{'L':>2} {'p_top':>8} {'p_bot':>8} {'dp_hPa':>7} {'DETA':>10}")
for L in range(50):
    print(f"{L+1:>2} {pint[L]:8.2f} {pint[L+1]:8.2f} {dp[L]:7.2f} {deta[L]:10.6f}")