// Sample-by-sample mechanical piano prototype. No recorded audio is loaded.
// Modal coordinates represent a stiff string with a moving bridge boundary.
// A felt hammer applies a unilateral nonlinear contact force at the strike point.
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <vector>
#include <algorithm>
#include <stdexcept>
constexpr double pi=3.14159265358979323846;
struct StringMode { double q=0,v=0,m,k,loss,a,strike,force=0; };
struct BoardMode { double q=0,v=0,m,k,loss,beta,radiation,acc=0; };
struct Geometry {double note,length,diameter,density;};
const Geometry geometry[]={ {21,2.05,.00110,.045},{36,1.45,.00105,.016},{48,1.03,.00099,.0084},{57,.80,.001,.006165},{60,.64,.00095,.005565},{72,.34,.00087,.004668},{84,.175,.00079,.003849},{96,.092,.00072,.003196},{108,.054,.00067,.002768} };
struct Voice {
 double sr,dt,hold,age=0,hx=0,hv,hm,stiffness,strikePosition,staticK=0,residualMass=0,denominator;
 int note,strings; bool hit=false; double peakForce=0,contactTime=0,initialEnergy=0;
 std::vector<StringMode> modes; std::vector<BoardMode> board;
 Voice(int n,double velocity,double duration,double rate,double feltScale=8,double bridgeMass=.16):sr(rate),dt(1/(rate*4)),hold(duration),note(n){
  int idx=0;while(idx<7&&n>geometry[idx+1].note)idx++;
  auto a=geometry[idx],b=geometry[idx+1];double u=std::clamp((n-a.note)/(b.note-a.note),0.,1.);
  auto blend=[&](double x,double y){return x*std::pow(y/x,u);};
  double L=blend(a.length,b.length),d=blend(a.diameter,b.diameter),mu=blend(a.density,b.density);
  double EI=2e11*pi*std::pow(d,4)/64;
  double f=440*std::pow(2.,(n-69)/12.);
  strings=n<33?1:n<57?2:3;
  hm=.010*std::pow(.30,std::clamp((n-21)/87.,0.,1.));hv=.25+3.8*std::pow(velocity,1.5);
  initialEnergy=.5*hm*hv*hv;
  stiffness=2e9*feltScale*std::pow(2.,(n-60)/36.);strikePosition=.115-.025*std::clamp((n-48)/48.,0.,1.);
  for(int string=0;string<strings;string++){
   double cents=strings==1?0:strings==2?(string==0?-.5:.5):(string-1)*.6;
   double freq=f*std::pow(2.,cents/1200.);
   double T=4*L*L*freq*freq*mu-EI*pi*pi/(L*L);
   if(T<=0)throw std::runtime_error("invalid string tension");
   staticK+=T/L;residualMass+=mu*L/3;
   for(int j=1;j<=512;j++){
    double omega=std::sqrt(T/mu*std::pow(j*pi/L,2)+EI/mu*std::pow(j*pi/L,4));
    if(omega/(2*pi)>14000)break;
    // Compensate the numerical oscillator's leapfrog dispersion.
    double numericalOmega=2/dt*std::sin(omega*dt/2);
    double m=mu*L/2,cross=mu*L*(j%2?1:-1)/(j*pi);
    modes.push_back({0,0,m,m*numericalOmega*numericalOmega,.12+.000008*std::pow(omega/(2*pi),1.5),cross,std::sin(j*pi*strikePosition)});
    residualMass-=cross*cross/m;
   }
  }
  // Orthotropic rectangular plate modes. Bracing raises effective bending
  // stiffness; dimensions, density and loss are model parameters, not an IR.
  double plateL=1.45,plateW=1.20,thickness=.009,rho=420,mass=rho*thickness*plateL*plateW/4;
  double Dx=1e10*std::pow(thickness,3)/12*8,Dy=.65e9*std::pow(thickness,3)/12*8,Dxy=std::sqrt(Dx*Dy)*.4;
  double bx=.28+.40*std::clamp((n-21)/87.,0.,1.),by=.39;
  for(int x=1;x<=8;x++)for(int y=1;y<=7;y++){
   double px=x*pi/plateL,py=y*pi/plateW;
   double omega=std::sqrt((Dx*std::pow(px,4)+2*Dxy*px*px*py*py+Dy*std::pow(py,4))/(rho*thickness));
   double hz=omega/(2*pi);if(hz>7000)continue;
   double beta=std::sin(x*pi*bx)*std::sin(y*pi*by),rad=std::sin(x*pi*.43)*std::sin(y*pi*.57);
   double Q=8+8*std::sqrt(hz/200.);
   board.push_back({0,0,mass,mass*omega*omega,omega/(2*Q),beta,rad});
  }
  residualMass+=bridgeMass;
  double inverseMass=0;for(auto &b:board)inverseMass+=b.beta*b.beta/b.m;
  denominator=1+residualMass*inverseMass;
 }
  double energy() const {
  double z=0,zv=0,E=.5*hm*hv*hv,fullMass=residualMass;
  for(auto &b:board){z+=b.beta*b.q;zv+=b.beta*b.v;E+=.5*b.m*b.v*b.v+.5*b.k*b.q*b.q;}
  double y=strikePosition*z;
  for(auto &m:modes){E+=.5*m.m*m.v*m.v+.5*m.k*m.q*m.q+m.a*m.v*zv;fullMass+=m.a*m.a/m.m;y+=m.strike*m.q/strings;}
  E+=.5*fullMass*zv*zv+.5*staticK*z*z;
  return E+stiffness*std::pow(std::max(0.,hx-y),3.3)/3.3;
 }
 double step(){
  double output=0;
  for(int sub=0;sub<4;sub++){
   double z=0,zv=0;for(auto &b:board){z+=b.beta*b.q;zv+=b.beta*b.v;}
   double contactY=strikePosition*z,contactV=strikePosition*zv;
   for(auto &m:modes){contactY+=m.strike*m.q/strings;contactV+=m.strike*m.v/strings;}
   double penetration=std::max(0.,hx-contactY);
   double F=stiffness*std::pow(penetration,2.3)*std::max(0.,1+.12*(hv-contactV));
   if(F>0){hit=true;contactTime+=dt;}peakForce=std::max(peakForce,F);
   if(!std::isfinite(F)||F>50000)throw std::runtime_error("unstable hammer contact");
   double damping=note>=89?0:age>hold?80*std::min(1.,(age-hold)/.025):0;
   double common=strikePosition*F-staticK*z;
   for(auto &m:modes){
    m.force=F/strings*m.strike-m.k*m.q-2*m.m*(m.loss+damping)*m.v;
    common-=m.a/m.m*m.force;
   }
   double betaAcc=0;
   for(auto &b:board){b.acc=(-b.k*b.q-2*b.m*b.loss*b.v+b.beta*common)/b.m;betaAcc+=b.beta*b.acc;}
   double correction=residualMass*betaAcc/denominator,zacc=0;
   for(auto &b:board){b.acc-=b.beta/b.m*correction;zacc+=b.beta*b.acc;}
   for(auto &m:modes){m.v+=dt*(m.force-m.a*zacc)/m.m;m.q+=dt*m.v;}
   hv-=dt*F/hm;hx+=dt*hv;
   double pressure=0;
   for(auto &b:board){b.v+=dt*b.acc;b.q+=dt*b.v;pressure+=b.radiation*b.v;}
   output+=pressure/4;age+=dt;
  }
  return output;
 }
};
int main(int argc,char**argv){
 if(argc<6||argc>8){std::fprintf(stderr,"usage: mechanics NOTE VELOCITY HOLD SECONDS output.f32 [felt-stiffness-scale] [bridge-mass-kg]\n");return 1;}
 try{
  int n=std::atoi(argv[1]);double velocity=std::atof(argv[2]),hold=std::atof(argv[3]),seconds=std::atof(argv[4]);
  Voice v(n,velocity,hold,44100,argc>=7?std::atof(argv[6]):8,argc==8?std::atof(argv[7]):.16);FILE*file=std::fopen(argv[5],"wb");if(!file)return 2;
  double peak=0;for(int i=0;i<int(seconds*44100);i++){float value=v.step();if(!std::isfinite(value))throw std::runtime_error("non-finite output");std::fwrite(&value,sizeof(value),1,file);peak=std::max(peak,std::abs(double(value)));}
  std::fclose(file);std::fprintf(stderr,"note=%d modes=%zu board=%zu peak=%g hammerHit=%d force=%g contactMs=%g energyRatio=%g\n",n,v.modes.size(),v.board.size(),peak,v.hit,v.peakForce,v.contactTime*1000,v.energy()/v.initialEnergy);
 }catch(const std::exception&e){std::fprintf(stderr,"%s\n",e.what());return 3;}
}