Skip to content

Instantly share code, notes, and snippets.

@kinoh kinoh/2mass.cpp
Last active Dec 19, 2015

Embed
What would you like to do?
IIshizaka and Flanagan's two-mass model of the vocal cords.
#include <iostream>
#include <cstdlib>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_odeiv.h>
#include <cmath>
const gsl_odeiv_step_type *StepType = gsl_odeiv_step_rk4;
const double dt = 1.E-5;
const double Hmin = 1.E-3;
const double T = 1.0;
const double m1 = 0.125;
const double m2 = 0.025;
const double ks1 = 80000;
const double ks2 = 8000;
const double es1 = 100;
const double es2 = 100;
const double kh1 = 3 * ks1;
const double kh2 = 3 * ks2;
const double eh1 = 500;
const double eh2 = 500;
const double kc = 25000;
const double om1 = std::sqrt(ks1 / m1);
const double om2 = std::sqrt(ks2 / m2);
const double xi1 = 0.1;
const double xi2 = 0.6;
const double lg = 1.4;
const double d1 = 0.25;
const double d2 = 0.05;
const double Ag0 = 0.05;
const double Ps = 7850;
const double xc = -Ag0 / (2 * lg);
//constexpr double A[] = { 0.45, 0.20, 0.26, 0.21, 0.32, 0.30, 0.33, 1.05, 1.12, 0.85, 0.63, 0.39, 0.26, 0.28, 0.23, 0.32, 0.29, 0.28, 0.40, 0.66, 1.20, 1.05, 1.62, 2.09, 2.56, 2.78, 2.86, 3.02, 3.75, 4.60, 5.09, 6.02, 6.55, 6.29, 6.27, 5.94, 5.28, 4.70, 3.87, 4.13, 4.25, 4.27, 4.69, 5.03 }; const char tract[] = "a";
//constexpr double A[] = { 0.33, 0.30, 0.36, 0.34, 0.68, 0.50, 2.43, 3.15, 2.66, 2.49, 3.39, 3.80, 3.78, 4.35, 4.50, 4.43, 4.68, 4.52, 4.15, 4.09, 3.51, 2.95, 2.03, 1.66, 1.38, 1.05, 0.60, 0.35, 0.32, 0.12, 0.10, 0.16, 0.25, 0.24, 0.38, 0.28, 0.36, 0.65, 1.58, 2.05, 2.01, 1.58 }; const char[] tract = "i";
//constexpr double A[] = { 0.40, 0.38, 0.28, 0.43, 0.55, 1.72, 2.91, 2.88, 2.37, 2.10, 3.63, 5.86, 5.63, 5.43, 4.80, 4.56, 4.29, 3.63, 3.37, 3.16, 3.31, 3.22, 2.33, 2.07, 2.07, 1.52, 0.74, 0.23, 0.15, 0.22, 0.22, 0.37, 0.60, 0.76, 0.86, 1.82, 2.35, 2.55, 3.73, 5.47, 4.46, 2.39, 1.10, 0.77, 0.41, 0.86 }; const char tract[] = "u";
//constexpr double A[] = { 0.21, 0.13, 0.16, 0.14, 0.06, 0.78, 1.25, 1.24, 0.99, 0.72, 0.73, 1.06, 1.77, 1.97, 2.46, 2.70, 2.92, 3.03, 2.84, 2.84, 2.83, 2.36, 2.14, 2.00, 1.78, 1.81, 1.79, 1.50, 1.37, 1.36, 1.43, 1.83, 2.08, 2.59, 2.54, 2.11, 2.34, 2.74, 2.19, 1.60 }; const char tract[] = "e";
constexpr double A[] = { 0.18, 0.17, 0.23, 0.28, 0.59, 1.46, 1.60, 1.11, 0.82, 1.01, 2.72, 2.71, 1.96, 1.92, 1.70, 1.66, 1.52, 1.28, 1.44, 1.28, 0.89, 1.25, 1.38, 1.09, 0.71, 0.46, 0.39, 0.32, 0.57, 1.06, 1.38, 2.29, 2.99, 3.74, 4.39, 5.38, 7.25, 7.00, 4.57, 2.75, 1.48, 0.68, 0.39, 0.14 }; const char tract[] = "o";
constexpr int TractN = sizeof(A) / sizeof(A[0]);
const double l = 0.397;
const double rho = 0.00114;
const double mu = 0.000186;
const double c = 35000;
constexpr double R(const int i) { return M_SQRTPI * std::sqrt(2. * rho * mu * om1) * std::pow(A[i], -1.5); }
constexpr double L(const int i) { return rho * l / (2. * A[i]); }
constexpr double C(const int i) { return l / (rho * c * c) * A[i]; }
constexpr double Ro(const int n) { return 128. * c * rho / (9. * M_PI * M_PI * A[n-1]); }
constexpr double Lo(const int n) { return 8. * rho / (3. * M_PI * std::sqrt(M_PI * A[n-1])); }
// y[] arrangement : [H1, V1, H2, V2, Ug, U1, ..., Un, P1, ..., Pn, Uo]
// int p[0] : current tract length
auto func = [](double t, const double y[], double dydt[], void *params) -> int
{
const int n = TractN;
const double h1 = Hmin > 0
? (y[0] > 0 && y[0] < Hmin ? Hmin : y[0])
: y[0];
const double v1 = y[1];
const double h2 = Hmin > 0
? (y[2] > 0 && y[2] < Hmin ? Hmin : y[2])
: y[2];
const double v2 = y[3];
const double Ug = y[4];
const double *const U = y + 5;
const double *const P = y + 5 + n;
const double Uo = y[5 + n * 2];
double *const dh1 = dydt + 0;
double *const dv1 = dydt + 1;
double *const dh2 = dydt + 2;
double *const dv2 = dydt + 3;
double *const dUg = dydt + 4;
double *const dU = dydt + 5;
double *const dP = dydt + 5 + n;
double *const dUo = dydt + 5 + n * 2;
const double x1 = h1 + xc;
const double x2 = h2 + xc;
const double s1 = ks1 * x1 * (1. + es1 * x1 * x1) + (h1 > 0 ? 0. : kh1 * h1 * (1. + eh1 * h1 * h1));
const double s2 = ks2 * x2 * (1. + es2 * x2 * x2) + (h2 > 0 ? 0. : kh2 * h2 * (1. + eh2 * h2 * h2));
const double r1 = 2. * m1 * om1 * (h1 > 0 ? xi1 : xi1 + 1.);
const double r2 = 2. * m2 * om2 * (h2 > 0 ? xi2 : xi2 + 1.);
const double Rbb = rho / (8. * lg * lg);
const double Rv = 1.5 * mu / lg;
const double Lg = rho / (2. * lg);
const double Rc = 1.37 * Rbb * Ug / (h1 * h1);
const double Rv1 = Rv * d1 / gsl_pow_3(h1);
const double Lg1 = Lg * d1 / h1;
const double R12 = Rbb * Ug * (1 / (h2 * h2) - 1 / (h1 * h1));
const double Rv2 = Rv * d2 / gsl_pow_3(h2);
const double Lg2 = Lg * d2 / h2;
const double Re = -0.5 * Rbb * Ug / (h2 * h2);
*dUg = (h1 <= 0 || h2 <= 0)
? -Ug / dt
: (Ps - P[0] - (Rc + Rv1 + R12 + Rv2 + Re) * Ug) / (Lg1 + Lg2);
// need dUg
const double dpi = -Rc * Ug;
const double dp1 = -Rv1 * Ug - Lg1 * *dUg;
const double dp12 = -R12 * Ug;
const double dp2 = -Rv2 * Ug - Lg2 * *dUg;
// const double dpo = Re * Ug;
const double p1 = Ps + (h1 <= 0 || h2 <= 0 ? 0. : dpi + dp1 / 2);
const double p2 = (h1 <= 0 ? 0. : Ps + (h2 <= 0 ? 0. : dpi + dp1 + dp12 + dp2 / 2));
for (int i = 0; i < n; i++)
dU[i] = (P[i] - P[i+1] - (R(i) + R(i+1)) * U[i]) / (L(i) + L(i+1));
dU[n-1] = (P[n-1] - R(n-1) * U[n-1]) / L(n-1);
*dUo = dU[n-1] - Ro(n) * Uo / Lo(n);
dP[0] = (Ug - U[0]) / C(0);
for (int i = 1; i < n; i++)
dP[i] = (U[i-1] - U[i]) / C(i);
*dh1 = v1;
*dv1 = (p1 * lg * d1 - kc * (h1 - h2) - s1 - r1 * v1) / m1;
*dh2 = v2;
*dv2 = (p2 * lg * d2 - kc * (h2 - h1) - s2 - r2 * v2) / m2;
return GSL_SUCCESS;
};
int main(int argc, char **argv)
{
if (argc - 1 < 5)
{
std::cerr << "Too Few Arguments;" << std::endl;
exit(EXIT_FAILURE);
}
const int d = 5 + TractN * 2 + 1;
int p[1];
gsl_odeiv_system sys = { func, NULL, d, p };
gsl_odeiv_step *step = gsl_odeiv_step_alloc(StepType, d);
double t = 0;
double y[d] = {0};
double dydt[d];
double yerr[d];
p[0] = sizeof(A) / sizeof(double);
y[0] = std::atof(argv[1]);
y[1] = std::atof(argv[2]);
y[2] = std::atof(argv[3]);
y[3] = std::atof(argv[4]);
y[4] = std::atof(argv[5]);
int m = (int)(T / dt);
for (int i = 0; i <= m; i++)
{
t = dt * i;
int status = gsl_odeiv_step_apply(step, t, dt, y, yerr, i == 0 ? NULL : dydt, dydt, &sys);
if (status)
{
std::cerr << "failed : " << status << std::endl;
exit(EXIT_FAILURE);
}
else if (gsl_isnan(yerr[0]))
{
std::cerr << "NaN error" << std::endl;
exit(EXIT_FAILURE);
}
std::printf("%.6f %11.4e %11.4e %11.4e %11.4e %11.4e %11.4e %11.4e %11.4e %11.4e\n",
t,
y[0],
y[1],
y[2],
y[3],
2. * lg * y[0],
2. * lg * y[2],
y[4],
y[5 + TractN * 2],
dydt[5 + TractN * 2]
);
}
return EXIT_SUCCESS;
}
@kinoh

This comment has been minimized.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.