31:T139b,animated = { const { L, h, initialOffset } = params; const { sourceOffset, sourceFreq } = srcParams; const { viewOptions } = viewParams; const { arrowScale, historyScale } = viewExpert; const xmax = Math.max(initialOffset, Math.abs(sourceOffset)); const ymax = h + L - Math.sqrt(L**2 - initialOffset**2) + 2; const Dt = arrowScale*dt, Dt2=0.5*Dt**2; // used to scale velocity vector const now = traj[when]; const x0 = now.x, y0 = now.y; // Compute fading trail of recent time steps const ntrail = 8 * historyScale; const trail = traj.slice(Math.max(0,when-ntrail), when+1); trail.forEach((t,k,T) => { const s = (k + 1) / (trail.length + 1); // [1/2] when length=1, [1/3,2/3] when length=2, etc T[k].alpha = s ** 3; }); const trail_nth = trail.filter((_,k) => k%historyScale == 0); const marks = [ // string Plot.link([now], { x1:"x", y1:"y", x2:0, y2:L+h, clip: true, stroke:"lightgray", strokeWidth:2 }), ]; if(viewOptions.includes("pos history")) marks.push( Plot.dot(trail_nth, { x:"x", y:"y", r:3, fill:"red", fillOpacity:"alpha" }) ); marks.push( // current position Plot.dot([now], { x:"x", y:"y", r:5, fill:"gray" }), ); if(viewOptions.includes("line of sight")) marks.push( Plot.link([now], { x1:"x", y1:"y", x2:sourceOffset, y2:0, stroke:"green", strokeDasharray:[6,2] }), ); if(viewOptions.includes("los velocity")) marks.push( Plot.arrow([now], { x1:"x", y1:"y", x2:d=>d.x+d.vxrel*Dt, y2:d=>d.y+d.vyrel*Dt, stroke:"green", strokeWidth:4, strokeOpacity:0.5, headLength: 4 }) ); if(viewOptions.includes("wavefronts")) { const wlen = speedOfSound / sourceFreq * 100; // in cm const oversampling = 24; const nfronts = Math.floor(Math.sqrt(ymax**2 + (sourceOffset+xmax)**2) / (wlen / oversampling)); /** const rdata = d3.range(1, nfronts+1).map(i => i * wlen / oversampling); console.log(rdata); marks.push(Plot.dot(rdata, { // ... channel definitions (x, y, etc.) x: r=>r, render: (index, scales, values) => { const xc = scales.x(sourceOffset), yc = scales.y(0); return htl.svg` ${Array.from(index, i => { // calculate arc parameters using scales and values[i] const r = scales.x(values.x[i]); const phase = (i % oversampling) / oversampling * 2 * Math.PI; const alpha = Math.cos(phase) ** 2; console.log(i, values.x[i], r, alpha); // ... other calculations //const pathData = `M${startX},${startY}A...`; // Construct SVG path data //return htl.svg``; return htl.svg``; })} `; }, }));**/ for(let i=1; i<=nfronts;i++) { const r = i * wlen / oversampling; const phase = (i % oversampling) / oversampling * 2 * Math.PI; const alpha = Math.cos(phase) ** 2; marks.push(Plot.line(halfCircle, { x:d=>sourceOffset+r*d.x, y:d=>r*d.y, clip: true, stroke:"green", strokeWidth:2, strokeOpacity:alpha, })); } } if(viewOptions.includes("vel history")) marks.push( Plot.link(trail_nth, { x1:x0, y1:y0, x2:d=>x0+d.vx*Dt, y2:d=>y0+d.vy*Dt, stroke:"blue", strokeOpacity:"alpha" }), ); if(viewOptions.includes("velocity")) marks.push( Plot.arrow([now], { x1:"x", y1:"y", x2:d=>d.x+d.vx*Dt, y2:d=>d.y+d.vy*Dt, stroke:"blue" }) ); if(viewOptions.includes("acceleration")) marks.push( Plot.arrow([now], { x1:"x", y1:"y", x2:d=>d.x+(d.aRx+d.aTx)*Dt2, y2:d=>d.y+(d.aRy+d.aTy)*Dt2, stroke:"magenta" }), ); if(viewOptions.includes("vel history") && viewOptions.includes("acceleration")) marks.push( Plot.line(trail, { x:d=>x0+d.vx*Dt, y:d=>y0+d.vy*Dt, stroke:"magenta", strokeOpacity:"alpha" }), Plot.arrow([now], { x1:d=>d.x+d.vx*Dt, y1:d=>d.y+d.vy*Dt, x2:d=>d.x+d.vx*Dt+(d.aRx+d.aTx)*Dt2, y2:d=>d.y+d.vy*Dt+(d.aRy+d.aTy)*Dt2, stroke:"magenta" }), ); if(viewOptions.includes("IMU accel")) marks.push( Plot.arrow([now], { x1:"x", y1:"y", x2:d=>d.x+(d.aXx+d.aYx)*Dt2, y2:d=>d.y+(d.aXy+d.aYy)*Dt2, stroke:"cyan" }), Plot.link([now], { x1:"x", y1:"y", x2:d=>d.x+d.aXx*Dt2, y2:d=>d.y+d.aXy*Dt2, stroke:"cyan" }), Plot.link([now], { x1:"x", y1:"y", x2:d=>d.x+d.aYx*Dt2, y2:d=>d.y+d.aYy*Dt2, stroke:"cyan" }) ); marks.push( // ground plane Plot.ruleY([0], { y:0, stroke:"black" }), // sound source Plot.dot([0], { x:sourceOffset, y:0, fill:"green", r:5 }), ); return Plot.plot({ marginTop: 8, marginRight: 8, x: { domain:[-xmax, +xmax], label: "x [cm]", labelAnchor: "center", labelArrow: false }, y: { domain:[0,ymax], label: "y [cm]", labelAnchor: "center", labelArrow: false }, aspectRatio: 1, width: width, grid: true, marks }); }32:T686,function makePlot(plot) { const marks = [ ]; let label = ""; if(plot) { const { name, unit, get, labels } = plot; label = `${name} [${unit}]`; for(let i=0; iget(d)[i], stroke:d=>labels[i], strokeWidth:4, strokeOpacity:0.5 }), Plot.dot([traj[when]], { x:"t", y:d=>get(d)[i], fill:d=>labels[i], r:4 }), ); } if(name === "IMU X-axis" && data !== null) marks.push( Plot.dot(data, { x:"t", y:d=>-d.aX, r:1.5, fill:"black", clip:true }), Plot.line(data, { x:"t", y:d=>-d.aX, stroke:"black", strokeOpacity:1, clip:true })); if(name === "Doppler Shift" && data !== null) marks.push( Plot.dot(data, { x:"t", y:d=>d.dfobs, r:1.5, fill:"black", clip:true }), Plot.line(data, { x:"t", y:d=>d.dfobs, stroke:"black", strokeOpacity:1, clip:true })); marks.push( // RenderFunction: https://github.com/observablehq/plot/blob/v0.6.17/src/mark.d.ts#L42-L60 (index, scales, values, dims, context) => { const width = 80; const height = 40; const x = dims.width - dims.marginRight - width; const y = dims.height - dims.marginBottom - height; return htl.svg`${Plot.legend(scales.scales)}`; }); } marks.push(Plot.frame()); return Plot.plot({ width: width, height: width/4, marginTop: 8, marginRight: 8, grid: true, x: { label: "Time [s]", labelAnchor: "center", labelArrow: false, domain:[0,tmax] }, y: { label, labelAnchor: "center", labelArrow: false }, color: { scheme:"Set1" }, marks }); }33:T902,function getTrajectory(params) { const { sqrt, sin, cos, PI } = Math; const { L, h, initialOffset, gamma } = params; const { A, dphase } = driveParams; const { sourceOffset, sourceFreq } = srcParams; const initialAngle = Math.asin(initialOffset / L); const phi = dphase * PI / 180; const omega0 = sqrt(g / L); params = Object.assign({g:981}, params); const yn = [[initialAngle,0]]; const steps = Math.ceil(tmax / dt); for (var i = 0; i < steps; i++) { yn.push(odeRK4([], yn[i], (yp, y, t) => { const theta = y[0]; const omega = y[1]; // d(theta)/dt = omega yp[0] = omega; // d(omega)/dt = -gamma * omega - (g / L) * sin(theta) yp[1] = -gamma * omega - (g / L) * sin(theta); // Add sinusoidal horizontal drive at pivot point yp[1] -= A * omega0**2 / L * sin(omega0 * t) * cos(theta); }, dt, i * dt)); } return yn.map(([theta,omega], i) => { const sth = sin(theta), cth = cos(theta); const rot2xy = (R,T) => [ R*sth + T*cth, -R*cth + T*sth ]; const omsq = omega**2; const x = L*sth, y = L + h - L*cth; const v = L*omega; const [ vx, vy ] = rot2xy(0,v); const srcdist = Math.sqrt((sourceOffset - x)**2 + y**2); const nx = (sourceOffset - x) / srcdist, ny = -y / srcdist; const vrel = vx * nx + vy * ny; const vxrel = vrel * nx, vyrel = vrel * ny; const dfobs = sourceFreq * (1 + vrel / (100*speedOfSound)) - sourceFreq; // radial and tangential components of coordinate acceleration in pivot inertial frame const theta_dotdot = -gamma*omega - (g/L)*sth; const aR = -L * omsq; const aT = L * theta_dotdot; // components of proper acceleration measured by IMU in its // rotating X (radial) and Y (tangential) frame const aX = L * omsq + g * cth; const aY = L * theta_dotdot + g * sth; // = -L * gamma * omega // convert from (radial,tangential) to (x,y) components const [aRx,aRy] = rot2xy(aR,0); const [aTx,aTy] = rot2xy(0,aT); const [aXx,aXy] = rot2xy(aX,0); const [aYx,aYy] = rot2xy(0,aY); ///const aX = g * cth + L * omsq, aY = -g * sth; return { t:dt*i, angle: theta * 180/PI, x, y, v, vx, vy, vrel, vxrel, vyrel, dfobs, aR, aT, aRx, aRy, aTx, aTy, aX, aY, aXx, aXy, aYx, aYy, }; }) }34:[] 2f:{"title":"Pendulum Doppler Shift","cells":[{"id":0,"value":"# Pendulum Doppler Shift","pinned":false,"hidden":false,"mode":"md"},{"id":129,"value":"viewof viewParams = Inputs.form({\n viewOptions: Inputs.checkbox([\n \"velocity\", \"acceleration\", \"IMU accel\", \"pos history\", \"vel history\",\n \"line of sight\", \"los velocity\", \"wavefronts\"\n ], {label:\"View Options\"}),\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":68,"value":"viewof when = Scrubber(d3.range(traj.length), { autoplay:false, loop:false, delay:null, format:d=>\"\" })","pinned":false,"hidden":false,"mode":"ojs"},{"id":48,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":482,"value":"viewof viewExpert = Inputs.form({\n arrowScale: Inputs.range([1, 20], { value:8, step:1, label:\"Arrow Scale\" }),\n historyScale: Inputs.range([1, 10], { value: 4, step:1, label:\"History Scale\" }),\n},\n{\n template: (inputs) => html`

${Object.values(inputs)}`\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":17,"value":"viewof params = Inputs.form({\n L: Inputs.range([50,200], { value:150, step:0.1, label:\"String Length [cm]\" }),\n h: Inputs.range([0,50], { value:10, step:0.1, label:\"Height [cm]\" }),\n initialOffset: Inputs.range([0,100], { value: 50, step:0.1, label:\"Initial Offset [cm]\" }),\n gamma: Inputs.range([0, .5], { value: 0.05, step:0.001, label:\"Damping Coef\" }),\n},\n{\n template: (inputs) => html`

${Object.values(inputs)}`\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":594,"value":"viewof driveParams = Inputs.form({\n A: Inputs.range([0,5], { value:0, step:0.01, label:\"Drive Amplitude [cm]\" }),\n dphase: Inputs.range([-180,180], { value:0, step:0.1, label:\"Drive Phase [deg]\" }),\n},\n{\n template: (inputs) => html`

${Object.values(inputs)}`\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":601,"value":"equilDrive = {\n const { L, initialOffset, gamma } = params;\n const initialAmp = Math.asin(initialOffset / L);\n const omega0 = Math.sqrt(g / L);\n return initialAmp * gamma * L / omega0;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":276,"value":"viewof plotSel = Inputs.form({\n plot1: Inputs.select(plots, { value:plots[7], label:\"Upper Plot\", format:plotLabel }),\n plot2: Inputs.select(plots, { value:plots[10], label:\"Lower Plot\", format:plotLabel }),\n},\n{\n template: (inputs) => html`

${Object.values(inputs)}`\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":35,"value":"plot1 = makePlot(plotSel.plot1)","pinned":false,"hidden":false,"mode":"ojs"},{"id":298,"value":"plot2 = makePlot(plotSel.plot2)","pinned":false,"hidden":false,"mode":"ojs"},{"id":188,"value":"viewof srcParams = Inputs.form({\n sourceOffset: Inputs.range([-100,100], { value: 25, step:0.1, label:\"Source Offset [cm]\" }),\n sourceFreq: Inputs.range([500,2500], { value: 2000, step:1, label:\"Source Freq [Hz]\" }),\n},\n{\n template: (inputs) => html`

${Object.values(inputs)}`\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":578,"value":"viewof soundParams = Inputs.form({\n airTemp: Inputs.range([0, 30], { value: 20, step:0.1, label:\"Air Temp [C]\" }),\n relHumidity: Inputs.range([0, 100], { value: 50, step:0.1, label:\"Rel Humidity [%]\" }),\n},\n{\n template: (inputs) => html`

${Object.values(inputs)}`\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":198,"value":"speedPrint = md`Speed of sound at ${soundParams.airTemp}°C and ${soundParams.relHumidity}% relative humidity is ${Number(speedOfSound).toFixed(2)} m/s`;","pinned":false,"hidden":false,"mode":"ojs"},{"id":367,"value":"viewof fileOptions = Inputs.form({\n file: Inputs.file({label: \"CSV file\", accept: \".csv\", required: true }),\n start: Inputs.range([0,1],{step:0.001, value:0.0, label:\"start\"}),\n stop: Inputs.range([0,1],{step:0.001, value:1.0, label:\"stop\"}),\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":378,"value":"dataPlot = {\n if(!allData) return Plot.plot();\n let tlo, thi;\n const tmax = allData[allData.length-1].t_ms;\n tlo = fileOptions.start * tmax;\n thi = fileOptions.stop * tmax;\n const what = d => -d.ax_g;\n return Plot.plot({\n width, \n height: width / 4,\n x: { label:\"Time [s]\", labelAnchor: \"center\", labelArrow: false },\n y: { label:\"IMU X Accel / g\", labelAnchor: \"center\", labelArrow: false },\n marginLeft: 64,\n marks: [\n Plot.line(allData?.filter(d => d.t_ms < tlo), { x:d=>d.t_ms/1000, y:what, stroke:\"gray\" }),\n Plot.line(allData?.filter(d => d.t_ms >= tlo && d.t_ms <= thi), {\n x:d=>d.t_ms/1000, y:what, stroke:\"blue\" }),\n Plot.line(allData?.filter(d => d.t_ms > thi), { x:d=>d.t_ms/1000, y:what, stroke:\"gray\" }),\n allData ? Plot.ruleX([tlo,thi], { x:t=>t/1000, stroke:\"red\" }) : null,\n Plot.frame(),\n ]\n });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":564,"value":"dataHistograms = html`\n

\n ${histPlot(d=>Math.sqrt(d.aX**2+d.aY**2), \"|IMU accel / g|\", \".4f\")}\n ${histPlot(d=>d.dfobs+srcParams.sourceFreq, \"Frequency [Hz]\", \".2f\")}\n

`","pinned":false,"hidden":false,"mode":"ojs"},{"id":337,"value":"md`---`","pinned":false,"hidden":false,"mode":"ojs"},{"id":517,"value":"function histPlot(what, label, prec) {\n if(data === null) return Plot.plot();\n const mean = d3.mean(data, what);\n const std = d3.deviation(data, what);\n const fmt = d3.format(prec);\n return Plot.plot({\n x: { label, labelAnchor: \"center\", labelArrow: false },\n marks: [\n Plot.rectY(data, Plot.binX({y: \"count\"}, {\n x: what, thresholds: 20,\n fill:\"blue\", fillOpacity:0.35, inset:0,\n })),\n Plot.ruleX([mean], { stroke:\"red\", strokeWidth:2 }),\n Plot.ruleX([mean-std,mean+std], { stroke:\"red\", strokeWidth:2, strokeDasharray:[3,6] }),\n Plot.text([`Mean = ${fmt(mean)}`], {\n fill:\"red\", x:mean, dx:-8, y:0, dy:-32, fontSize:20, textAnchor:\"end\"\n }),\n Plot.text([`Std.Dev. = ${fmt(std)}`], {\n fill:\"red\", x:mean+std, dx:8, y:0, dy:-64, fontSize:20, textAnchor:\"start\"\n }),\n Plot.frame(),\n ]\n });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":421,"value":"data = {\n if(!allData) return null;\n const tmax = allData[allData.length-1].t_ms;\n const tlo = fileOptions.start * tmax;\n const thi = fileOptions.stop * tmax;\n return allData.filter(d => d.t_ms >= tlo && d.t_ms <= thi).map(d => ({\n t: (d.t_ms - tlo)/1000,\n aX: d.ax_g,\n aY: Math.sqrt(d.ay_g ** 2 + d.az_g ** 2),\n dfobs: d.freq_hz - srcParams.sourceFreq,\n }));\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":374,"value":"allData = fileOptions.file?.csv({typed:true}) || null","pinned":false,"hidden":false,"mode":"ojs"},{"id":467,"value":"plotLabel = d => d === null ? \"\" : `${d.name} [${d.unit}]`","pinned":false,"hidden":false,"mode":"ojs"},{"id":278,"value":"plots = [\n null,\n { name: \"Angle\", unit: \"deg\", get:d=>[d.angle], labels:[\"θ\"] },\n { name: \"Position\", unit: \"cm\", get:d=>[d.x,d.y], labels:[\"x\",\"y\"] },\n { name: \"Velocity\", unit: \"cm/s\", get:d=>[d.vx,d.vy], labels:[\"vx\",\"vy\"] },\n { name: \"Coord Accel (x,y)\", unit: \"cm/s^2\", get:d=>[d.aRx+d.aTx,d.aRy+d.aTy], labels:[\"ax\",\"ay\"] },\n { name: \"Coord Accel (R,T)\", unit: \"cm/s^2\", get:d=>[d.aR,d.aT], labels:[\"aR\",\"aT\"] },\n { name: \"Proper Accel\", unit: \"cm/s^2\", get:d=>[d.aXx+d.aYx,d.aXy+d.aYy+g], labels:[\"Ax\",\"Ay+g\"] },\n { name: \"IMU X-axis\", unit: \"g\", get:d=>[d.aX/g], labels:[\"IMU-X\"] },\n { name: \"IMU Y-axis\", unit: \"g\", get:d=>[d.aY/g], labels:[\"IMU-Y\"] },\n { name: \"LOS Speed\", unit: \"cm/s\", get:d=>[d.vrel], labels:[\"vlos\"] },\n { name: \"Doppler Shift\", unit: \"Hz\", get:d=>[d.dfobs], labels:[\"Δf\"] },\n]","pinned":false,"hidden":false,"mode":"ojs"},{"id":286,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":194,"value":"speedOfSound = {\n const { airTemp, relHumidity } = soundParams;\n const P_atm = 1013.25; // Standard pressure in hPa\n \n // 1. Calculate speed of sound in dry air\n const c_dry = 331.3 * Math.sqrt(1 + airTemp / 273.15);\n \n // 2. Calculate partial pressure of water vapor (Saturation pressure * RH)\n // Magnus-Tetens Formula for Saturation Vapor Pressure\n const p_sat = 6.112 * Math.exp((17.62 * airTemp) / (243.12 + airTemp));\n const p_w = (relHumidity / 100) * p_sat;\n \n // 3. Apply humidity correction\n // The ratio of specific heats and gas constants leads to this multiplier\n const humidityMultiplier = Math.sqrt(1 + 0.511 * (p_w / P_atm));\n \n return c_dry * humidityMultiplier;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":174,"value":"halfCircle = {\n const { sin, cos, PI } = Math;\n const npts = 64;\n const dphi = PI / (npts - 1);\n return d3.range(npts).map(k => {\n const phi = k * dphi;\n return { x: cos(phi), y: sin(phi) };\n });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":33,"value":"traj = getTrajectory(params)","pinned":false,"hidden":false,"mode":"ojs"},{"id":5,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":319,"value":"g = 981","pinned":false,"hidden":false,"mode":"ojs"},{"id":52,"value":"dt = 1 / 60","pinned":false,"hidden":false,"mode":"ojs"},{"id":50,"value":"tmax = 30","pinned":false,"hidden":false,"mode":"ojs"},{"id":64,"value":"import {Scrubber} from \"@mbostock/scrubber\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":8,"value":"import {odeRK4} from \"@rreusser/integration\"","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W34","stdlib":"1","titling":"auto"}