"""Transient absolute-acceleration SRS reference. SI internally; no cargo release decision."""
import math,json,argparse
from pathlib import Path
G=9.80665

def spectrum(samples,fs,frequencies,zeta=.05,step_limit=.15):
 if not isinstance(samples,list) or len(samples)<2 or len(samples)>200000:raise ValueError('Supply 2..200000 acceleration samples in g')
 if any(isinstance(v,bool) or not isinstance(v,(int,float)) or not math.isfinite(v) for v in samples):raise ValueError('Finite numeric acceleration samples required')
 if not math.isfinite(fs) or fs<=0 or not .01<=zeta<=1 or not .02<=step_limit<=.2:raise ValueError('Invalid sampling, damping or integration setting')
 if not frequencies or len(frequencies)>100 or any(not math.isfinite(f) or f<=0 or f>fs/10 for f in frequencies):raise ValueError('Positive frequencies at most fs/10 required')
 output=[]
 for f in frequencies:
  w=2*math.pi*f;n=max(1,math.ceil(w/(fs*step_limit)));h=1/(fs*n)
  tail=math.ceil(5/(zeta*w)*fs);steps=(len(samples)-1+tail)*n
  if steps>2000000:raise ValueError('Reference resource limit exceeded')
  z=v=positive=negative=0.
  def deriv(x,y,a):return y,-2*zeta*w*y-w*w*x-a
  for i in range(len(samples)-1+tail):
   a0=samples[i]*G if i<len(samples) else 0.
   a1=samples[i+1]*G if i+1<len(samples) else 0.
   for k in range(n):
    a=a0+(a1-a0)*k/n;b=a0+(a1-a0)*(k+.5)/n;c=a0+(a1-a0)*(k+1)/n
    k1=deriv(z,v,a);k2=deriv(z+h*k1[0]/2,v+h*k1[1]/2,b);k3=deriv(z+h*k2[0]/2,v+h*k2[1]/2,b);k4=deriv(z+h*k3[0],v+h*k3[1],c)
    z+=h*(k1[0]+2*k2[0]+2*k3[0]+k4[0])/6;v+=h*(k1[1]+2*k2[1]+2*k3[1]+k4[1])/6
    response=(-2*zeta*w*v-w*w*z)/G
    positive=max(positive,response);negative=min(negative,response)
  output.append({'natural_frequency_hz':f,'positive_g':positive,'negative_g':negative,'absolute_g':max(positive,-negative),'ringdown_seconds':tail/fs,'integration_steps':steps})
 return {'response_quantity':'absolute acceleration','damping_ratio':zeta,'quality_factor':1/(2*zeta),'sampling_hz':fs,'input_unit':'g, standard gravity 9.80665 m/s^2','input_peak_g':max(abs(x) for x in samples),'results':output,'scope':'Linear SDOF, linearly interpolated measured input, five decay-time-constant ringdown; not a shock qualification or an anti-alias filter'}

def pulse(fs=1200,seconds=.01,peak=1):
 return [peak*math.sin(math.pi*i/fs/seconds) if i/fs<=seconds else 0. for i in range(math.ceil(seconds*fs)+2)]

def main():
 p=argparse.ArgumentParser();p.add_argument('--input');p.add_argument('--output');a=p.parse_args()
 data=json.loads(Path(a.input).read_text(encoding='utf-8-sig')) if a.input else {'samples':pulse(),'fs':1200,'frequencies':[20,40,80]}
 r=spectrum(**data)
 if a.output:Path(a.output).write_text(json.dumps(r,indent=2),encoding='utf-8')
 print(json.dumps(r,indent=2))
if __name__=='__main__':main()