oscillator.h 2.47 KB
 Dmitriy Morozov committed Feb 14, 2016 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 ``````#pragma once #include #include #include #include #include struct Oscillator { using Vertex = grid::Point; static constexpr float pi = 3.14159265358979323846; inline float evaluate(Vertex v, float t) const { t *= 2*pi; float dist2 = (center - v).norm(); float dist_damp = exp(-dist2/(2*radius*radius)); if (type == damped) { float phi = acos(zeta); float val = 1. - exp(-zeta*omega0*t) * (sin(sqrt(1-zeta*zeta)*omega0*t + phi) / sin(phi)); return val * dist_damp; } else if (type == decaying) { t += 1. / omega0; float val = sin(t / omega0) / (omega0 * t); return val * dist_damp; } else if (type == periodic) { t += 1. / omega0; float val = sin(t / omega0); return val * dist_damp; } else return 0; // impossible } Vertex center; float radius; float omega0; float zeta; enum { damped, decaying, periodic } type; }; using Oscillators = std::vector; // trim() from http://stackoverflow.com/a/217605/44738 static inline std::string <rim(std::string &s) { s.erase(s.begin(), std::find_if(s.begin(), s.end(), std::not1(std::ptr_fun(std::isspace)))); return s; } static inline std::string &rtrim(std::string &s) { s.erase(std::find_if(s.rbegin(), s.rend(), std::not1(std::ptr_fun(std::isspace))).base(), s.end()); return s; } static inline std::string &trim(std::string &s) { return ltrim(rtrim(s)); } Oscillators read_oscillators(std::string fn) { Oscillators res; std::ifstream in(fn); if (!in) throw std::runtime_error("Unable to open " + fn); std::string line; while(std::getline(in, line)) { line = trim(line); if (line.empty() || line[0] == '#') continue; std::istringstream iss(line); std::string stype; iss >> stype; auto type = Oscillator::periodic; if (stype == "damped") type = Oscillator::damped; else if (stype == "decaying") type = Oscillator::decaying; int x,y,z; iss >> x >> y >> z; float r, omega0, zeta; iss >> r >> omega0; if (type == Oscillator::damped) iss >> zeta; res.emplace_back(Oscillator { {x,y,z}, r, omega0, zeta, type }); } return res; }``````