import math import numpy as np from numba import njit from probes import OBSERVER, simulate as frozen_simulate from complex_bank import solve_bank energy = OBSERVER.ENGINE.energy accumulate = OBSERVER.ENGINE.accumulate smoothstep = OBSERVER.ENGINE.smoothstep @njit(cache=True) def integrate(masses, springs, pick_weights, pickup_weights, bridge_mass, bridge_stiffness, bridge_damping, port_weights, fret_nodes, event_start, event_time, event_string, event_pitch, event_velocity, event_picked, event_gate, event_contacts, pick_order, open_pitches, internal_loss, mass_loss, rate, seconds, finger_stiffness, finger_ratio, hammer_seconds, pull_seconds, prepare_seconds, finger_travel, pick_stiffness, pick_damping, pick_travel, hook_deflection, stroke_seconds, release_seconds, mute_rate, loss_scale, exact_release_time=True, consistent_motion=True): strings, nodes = masses.shape modes = len(port_weights) frames = math.ceil(seconds * rate) time_step = 1 / rate half = time_step / 2 root = math.sqrt(3) shift = time_step * (.25 + 1j / (4 * root)) square = shift ** 2 position = np.zeros((strings, nodes)) speed = np.zeros((strings, nodes)) bridge_position = np.zeros(modes) bridge_speed = np.zeros(modes) diagonal = np.empty((strings, nodes), dtype=np.complex128) off_diagonal = np.empty((strings, nodes - 1), dtype=np.complex128) rhs = np.empty((strings, nodes), dtype=np.complex128) bridge_rhs = np.empty(modes, dtype=np.complex128) bridge_diagonal = bridge_mass + shift * bridge_damping + square * bridge_stiffness coupling = np.empty(strings, dtype=np.complex128) pivots = np.empty((strings, nodes), dtype=np.complex128) reduced = np.empty((strings, nodes), dtype=np.complex128) free_solution = np.empty((strings, nodes), dtype=np.complex128) boundary_response = np.empty((strings, nodes), dtype=np.complex128) pick_response = np.empty((strings, nodes), dtype=np.complex128) acceleration = np.empty((strings, nodes), dtype=np.complex128) stage_position = np.empty((2, strings, nodes)) stage_velocity = np.empty((2, strings, nodes)) stage_bridge_position = np.empty((2, modes)) stage_bridge_velocity = np.empty((2, modes)) bridge_acceleration = np.empty(modes, dtype=np.complex128) fractions = np.zeros((strings, 24)) prior_fractions = np.zeros((strings, 24)) destination_fractions = np.zeros((strings, 24)) finger_target = np.zeros((strings, 24)) finger_velocity = np.zeros((strings, 24)) finger_damping = np.empty((strings, 24)) current_pitch = open_pitches.copy() old_pitch = open_pitches.copy() finger_begin = np.full(strings, -1.) finger_kind = np.zeros(strings, dtype=np.int64) finger_distance = np.zeros(strings) gate = np.zeros(strings) losses = np.empty(strings) pick_begin = np.full(strings, -1.) pick_release = np.full(strings, -1.) pick_origin = np.zeros(strings) pick_distance = np.zeros(strings) pick_fraction = np.zeros(strings) previous_pick_stiffness = np.zeros(strings) pick_target = np.zeros(strings) pick_velocity = np.zeros(strings) pick_coefficient = np.zeros(strings, dtype=np.complex128) output = np.empty((frames, strings)) stride = max(1, rate // 100) ledger = np.empty((2 + frames // stride, 11)) ledger[0] = np.zeros(11) samples = 1 finger_work = 0. pick_work = 0. dissipated = 0. finger_dissipated = 0. pick_dissipated = 0. bridge_exchange = 0. discharged_pick_energy = 0. maximum_balance_error = 0. maximum_pick_force = 0. compensation = np.zeros(7) pick_log = np.empty((len(event_time), 9)) pick_indices = np.full(strings, -1, dtype=np.int64) picks_logged = 0 finger_event = 0 pick_event = 0 for string in range(strings): coupling[string] = springs[string, -1] * (square + shift * internal_loss[string]) for fret in range(24): finger_damping[string, fret] = 2 * finger_ratio * math.sqrt(finger_stiffness * masses[string, fret_nodes[fret + 1]]) for node in range(nodes - 1): off_diagonal[string, node] = -(square + shift * internal_loss[string]) * springs[string, node + 1] for frame in range(frames): time = (frame + .5) / rate step_finger_work = 0. step_pick_work = 0. step_dissipated = 0. step_finger_dissipated = 0. step_pick_dissipated = 0. step_bridge_exchange = 0. step_discharged = 0. while finger_event < len(event_start) and time >= event_start[finger_event]: string = event_string[finger_event] prior_fractions[string] = fractions[string] destination_fractions[string] = event_contacts[finger_event] old_pitch[string] = current_pitch[string] current_pitch[string] = event_pitch[finger_event] finger_begin[string] = event_start[finger_event] finger_kind[string] = 0 if event_picked[finger_event] else 1 if current_pitch[string] > old_pitch[string] else -1 finger_distance[string] = finger_travel * event_velocity[finger_event] gate[string] = event_gate[finger_event] finger_event += 1 for string in range(strings): duration = prepare_seconds if finger_kind[string] == 0 else hammer_seconds if finger_kind[string] == 1 else pull_seconds phase = (time - finger_begin[string]) / duration mix = smoothstep(phase) for fret in range(24): fractions[string, fret] = prior_fractions[string, fret] * (1 - mix) + destination_fractions[string, fret] * mix finger_target[string, fret] = 0. finger_velocity[string, fret] = 0. hammer_active = finger_kind[string] == 1 and consistent_motion and (time - half - finger_begin[string]) / duration < 1 if finger_begin[string] >= 0 and (phase < 1 or hammer_active) and finger_kind[string] != 0: contact = current_pitch[string] - open_pitches[string] - 1 if finger_kind[string] == 1 else old_pitch[string] - open_pitches[string] - 1 if contact >= 0: if finger_kind[string] == 1: finger_target[string, contact] = finger_distance[string] * (1 - phase) finger_velocity[string, contact] = -finger_distance[string] / duration if consistent_motion: start_phase = min(1., max(0., (time - half - finger_begin[string]) / duration)) end_phase = min(1., max(0., (time + half - finger_begin[string]) / duration)) start_target = finger_distance[string] * (1 - start_phase) end_target = finger_distance[string] * (1 - end_phase) finger_target[string, contact] = .5 * (start_target + end_target) finger_velocity[string, contact] = (end_target - start_target) / time_step else: movement = min(1., phase / .6) finger_target[string, contact] = -finger_distance[string] * smoothstep(movement) if movement < 1: finger_velocity[string, contact] = -finger_distance[string] * 6 * movement * (1 - movement) / (.6 * duration) fractions[string, contact] = prior_fractions[string, contact] * (1 - smoothstep((phase - .6) / .4)) mute = min(1., max(0., (time - gate[string] + .024) / .024)) losses[string] = mass_loss + mute_rate * loss_scale * mute * mute while pick_event < len(event_time) and time >= event_time[pick_order[pick_event]]: selected_event = pick_order[pick_event] if event_picked[selected_event]: string = event_string[selected_event] pick_begin[string] = event_time[selected_event] pick_release[string] = -1. pick_origin[string] = 0. velocity_at_pick = 0. for node in range(nodes): pick_origin[string] += pick_weights[string, node] * position[string, node] velocity_at_pick += pick_weights[string, node] * speed[string, node] direction = 1. if picks_logged % 2 == 0 else -1. pick_distance[string] = direction * pick_travel * event_velocity[selected_event] pick_indices[string] = picks_logged pick_log[picks_logged] = np.array([time, float(string), direction, pick_origin[string], velocity_at_pick, pick_distance[string], energy(position, speed, bridge_position, bridge_speed, masses, springs, bridge_mass, bridge_stiffness, port_weights), -1., 0.]) picks_logged += 1 pick_event += 1 bridge_port = 0. bridge_velocity = 0. for mode in range(modes): bridge_port += port_weights[mode] * bridge_position[mode] bridge_velocity += port_weights[mode] * bridge_speed[mode] for string in range(strings): phase = (time - pick_begin[string]) / stroke_seconds pick_target[string] = pick_origin[string] + pick_distance[string] * smoothstep(phase) pick_velocity[string] = 0. if 0 <= phase < 1: pick_velocity[string] = pick_distance[string] * 6 * phase * (1 - phase) / stroke_seconds if consistent_motion and pick_begin[string] >= 0: start_phase = (time - half - pick_begin[string]) / stroke_seconds end_phase = (time + half - pick_begin[string]) / stroke_seconds start_target = pick_origin[string] + pick_distance[string] * smoothstep(start_phase) end_target = pick_origin[string] + pick_distance[string] * smoothstep(end_phase) pick_target[string] = .5 * (start_target + end_target) pick_velocity[string] = (end_target - start_target) / time_step pick_position = 0. for node in range(nodes): pick_position += pick_weights[string, node] * position[string, node] if pick_begin[string] >= 0 and pick_release[string] < 0: if abs(pick_target[string] - pick_position) >= hook_deflection or phase >= 1: pick_release[string] = pick_begin[string] + stroke_seconds if exact_release_time and phase >= 1 else time pick_log[pick_indices[string], 7] = pick_release[string] pick_log[pick_indices[string], 8] = 2. if phase >= 1 else 1. fraction = 0. if pick_begin[string] >= 0: fraction = 1. if pick_release[string] >= 0: release_phase = min(1., (time - pick_release[string]) / release_seconds) fraction = .5 * (1 + math.cos(math.pi * release_phase)) pick_fraction[string] = fraction stiffness = pick_stiffness * fraction damping = pick_damping * fraction pick_coefficient[string] = square * stiffness + shift * damping if stiffness < previous_pick_stiffness[string]: step_discharged += .5 * (previous_pick_stiffness[string] - stiffness) * (pick_position - pick_target[string]) ** 2 previous_pick_stiffness[string] = stiffness for node in range(nodes): left_position = position[string, node - 1] if node else 0. right_position = position[string, node + 1] if node + 1 < nodes else bridge_port left_velocity = speed[string, node - 1] if node else 0. right_velocity = speed[string, node + 1] if node + 1 < nodes else bridge_velocity spring_force = springs[string, node] * (position[string, node] - left_position) spring_force += springs[string, node + 1] * (position[string, node] - right_position) stiffness_speed = springs[string, node] * (speed[string, node] - left_velocity) stiffness_speed += springs[string, node + 1] * (speed[string, node] - right_velocity) viscous = 2 * losses[string] * masses[string, node] * speed[string, node] viscous += internal_loss[string] * stiffness_speed rhs[string, node] = -viscous - spring_force - shift * stiffness_speed diagonal[string, node] = masses[string, node] * (1 + 2 * shift * losses[string]) diagonal[string, node] += (square + shift * internal_loss[string]) * (springs[string, node] + springs[string, node + 1]) for fret in range(24): node = fret_nodes[fret + 1] stiffness = finger_stiffness * fractions[string, fret] damping = finger_damping[string, fret] * fractions[string, fret] diagonal[string, node] += square * stiffness + shift * damping rhs[string, node] += stiffness * (finger_target[string, fret] - position[string, node]) rhs[string, node] += damping * (finger_velocity[string, fret] - speed[string, node]) rhs[string, node] -= shift * stiffness * speed[string, node] pick_speed = 0. for node in range(nodes): pick_speed += pick_weights[string, node] * speed[string, node] forcing = pick_stiffness * fraction * (pick_target[string] - pick_position) forcing += pick_damping * fraction * (pick_velocity[string] - pick_speed) forcing -= shift * pick_stiffness * fraction * pick_speed for node in range(nodes): rhs[string, node] += pick_weights[string, node] * forcing for mode in range(modes): spring_force = bridge_stiffness[mode] * bridge_position[mode] stiffness_speed = bridge_stiffness[mode] * bridge_speed[mode] viscous = bridge_damping[mode] * bridge_speed[mode] for string in range(strings): spring_force += port_weights[mode] * springs[string, -1] * (bridge_port - position[string, -1]) link_speed = port_weights[mode] * springs[string, -1] * (bridge_velocity - speed[string, -1]) stiffness_speed += link_speed viscous += internal_loss[string] * link_speed bridge_rhs[mode] = -viscous - spring_force - shift * stiffness_speed solve_bank(diagonal, off_diagonal, bridge_diagonal, port_weights, coupling, pick_coefficient, pick_weights, rhs, bridge_rhs, pivots, reduced, free_solution, boundary_response, pick_response, acceleration, bridge_acceleration) for string in range(strings): for node in range(nodes): derivative = speed[string, node] + shift * acceleration[string, node] first_velocity = derivative.real - (2 - root) * derivative.imag second_velocity = derivative.real + (2 + root) * derivative.imag stage_velocity[0, string, node] = first_velocity stage_velocity[1, string, node] = second_velocity stage_position[0, string, node] = position[string, node] + time_step * (.25 * first_velocity + (.25 - root / 6) * second_velocity) stage_position[1, string, node] = position[string, node] + time_step * ((.25 + root / 6) * first_velocity + .25 * second_velocity) for mode in range(modes): derivative = bridge_speed[mode] + shift * bridge_acceleration[mode] first_velocity = derivative.real - (2 - root) * derivative.imag second_velocity = derivative.real + (2 + root) * derivative.imag stage_bridge_velocity[0, mode] = first_velocity stage_bridge_velocity[1, mode] = second_velocity stage_bridge_position[0, mode] = bridge_position[mode] + time_step * (.25 * first_velocity + (.25 - root / 6) * second_velocity) stage_bridge_position[1, mode] = bridge_position[mode] + time_step * ((.25 + root / 6) * first_velocity + .25 * second_velocity) for stage in range(2): bridge_velocity_port = sum(port_weights * stage_bridge_velocity[stage]) bridge_position_port = sum(port_weights * stage_bridge_position[stage]) for mode in range(modes): step_dissipated += half * bridge_damping[mode] * stage_bridge_velocity[stage, mode] ** 2 for string in range(strings): for fret in range(24): node = fret_nodes[fret + 1] velocity = stage_velocity[stage, string, node] contact_position = stage_position[stage, string, node] stiffness = finger_stiffness * fractions[string, fret] damping = finger_damping[string, fret] * fractions[string, fret] force = stiffness * (finger_target[string, fret] - contact_position) force += damping * (finger_velocity[string, fret] - velocity) step_finger_work += half * force * velocity step_finger_dissipated += half * damping * (velocity - finger_velocity[string, fret]) ** 2 pick_position_port = 0. pick_velocity_port = 0. for node in range(nodes): pick_position_port += pick_weights[string, node] * stage_position[stage, string, node] pick_velocity_port += pick_weights[string, node] * stage_velocity[stage, string, node] stiffness = pick_stiffness * pick_fraction[string] damping = pick_damping * pick_fraction[string] force = stiffness * (pick_target[string] - pick_position_port) force += damping * (pick_velocity[string] - pick_velocity_port) maximum_pick_force = max(maximum_pick_force, abs(force)) step_pick_work += half * force * pick_velocity_port step_pick_dissipated += half * damping * (pick_velocity_port - pick_velocity[string]) ** 2 previous_velocity = 0. for node in range(nodes): velocity = stage_velocity[stage, string, node] step_dissipated += half * 2 * losses[string] * masses[string, node] * velocity ** 2 step_dissipated += half * internal_loss[string] * springs[string, node] * (velocity - previous_velocity) ** 2 previous_velocity = velocity step_dissipated += half * internal_loss[string] * springs[string, -1] * (bridge_velocity_port - previous_velocity) ** 2 bridge_force = springs[string, -1] * (stage_position[stage, string, -1] - bridge_position_port) bridge_force += internal_loss[string] * springs[string, -1] * (previous_velocity - bridge_velocity_port) step_bridge_exchange += half * abs(bridge_force * bridge_velocity_port) finger_work, compensation[0] = accumulate(finger_work, compensation[0], step_finger_work) pick_work, compensation[1] = accumulate(pick_work, compensation[1], step_pick_work) dissipated, compensation[2] = accumulate(dissipated, compensation[2], step_dissipated) finger_dissipated, compensation[3] = accumulate(finger_dissipated, compensation[3], step_finger_dissipated) pick_dissipated, compensation[4] = accumulate(pick_dissipated, compensation[4], step_pick_dissipated) bridge_exchange, compensation[5] = accumulate(bridge_exchange, compensation[5], step_bridge_exchange) discharged_pick_energy, compensation[6] = accumulate(discharged_pick_energy, compensation[6], step_discharged) for string in range(strings): for node in range(nodes): position[string, node] += half * (stage_velocity[0, string, node] + stage_velocity[1, string, node]) speed[string, node] += time_step * (acceleration[string, node].real + root * acceleration[string, node].imag) for mode in range(modes): bridge_position[mode] += half * (stage_bridge_velocity[0, mode] + stage_bridge_velocity[1, mode]) bridge_speed[mode] += time_step * (bridge_acceleration[mode].real + root * bridge_acceleration[mode].imag) bridge_readout = sum(port_weights * bridge_speed) for string in range(strings): signal = -bridge_readout for node in range(nodes): signal += pickup_weights[string, node] * speed[string, node] output[frame, string] = signal if frame % stride == stride - 1 or frame == frames - 1: stored = energy(position, speed, bridge_position, bridge_speed, masses, springs, bridge_mass, bridge_stiffness, port_weights) contact_storage = 0. for string in range(strings): pick_position = 0. for node in range(nodes): pick_position += pick_weights[string, node] * position[string, node] contact_storage += .5 * pick_stiffness * pick_fraction[string] * (pick_position - pick_target[string]) ** 2 for fret in range(24): contact_storage += .5 * finger_stiffness * fractions[string, fret] * (position[string, fret_nodes[fret + 1]] - finger_target[string, fret]) ** 2 residual = stored - finger_work - pick_work + dissipated maximum_balance_error = max(maximum_balance_error, abs(residual)) ledger[samples] = np.array([(frame + 1) / rate, stored, stored + contact_storage, finger_work, pick_work, dissipated, finger_dissipated, pick_dissipated, residual, bridge_exchange, discharged_pick_energy]) samples += 1 return output, ledger[:samples], pick_log[:picks_logged], maximum_balance_error, maximum_pick_force, position, speed, bridge_position, bridge_speed def simulate(events, seconds, release=.00015, *, exact_release_time=True, consistent_motion=True, **controls): previous = OBSERVER.ENGINE.integrate try: OBSERVER.ENGINE.integrate = lambda *arguments: integrate(*arguments, exact_release_time, consistent_motion) return frozen_simulate(events, seconds, release, **controls) finally: OBSERVER.ENGINE.integrate = previous