have something happen every (e.g.) 15th cycle of an oscillator

[bgcm]

Write your question here.
here is code that works every cycle:
// non-lin (oscillator), Last point approx. (Cromer) version (w/ quadratic dissipation and momentary drive).
// Is an approximation for the numerical integration of a clock's equation of motion: ~ a Synchronome

#include <iostream>
#include <fstream> /* for output file */
#include <cmath>
#include <iomanip> /* for setw() */
#include <stdlib.h> /* for exit() */
using namespace std;

double const pi = 4*atan(1);

////////////
// printout routine
//
void p(ofstream& ouch,
double const t,
double const theta,
double const thetadot,
double const deltag, // ((4*pi*pi)/(period*period))*(1+deltag_now)
double const dr, // drive
double const r, // dissipation
double const energy

){

ouch << setprecision(10) << setw(15) << t << "," //the "," is CSV. Can be tabs?

//ouch << setw(10) << t << ","
<< setw(15) << theta << ","
<< setw(15) << thetadot << ","
<< setw(15) << deltag << ","
<< setw(15) << dr << ","
<< setw(15) << r << ","
<< setw(15) << energy << ","
<< endl;
}

// The expression "\"" is just too ugly.
// Define a nicer-looking synonym:
string const Q = "\"";

// C++ modulo operator "%" only works for integers
// Generalize to doubles
double mod(double const num, double const denom){
return num - denom * floor(num / denom);
}

////////////////////////////////////////
// Main program
//
int main() {
cout << "/////////////////////////////////////" << endl;
cout << "Is the numerical model of a clock pendulum," << endl;
cout << "w/ gravity type drive every 30s, and quadratic dissipation. " << endl;
cout << "and uses the LPA aproximation." << endl;
cout << "///////////////////////////////" << endl;

double period = 2;
double theta = 0.1; //radian ~ 4" / one meter
double d = 1e-4; //temporarily
double dr = 0;
double thetadot = 0;
double energy = thetadot*thetadot * 100;
double t = 0;
double steps_per_cycle = 1e3;
double deltat = period / steps_per_cycle;
double r =3.6e-4; // temporary 'till driving set

cout << "r ~= 3.6e-4 means Q~= 8.7k for a 2s period" << endl;
// cout << "enter drive amplitude; try 1.0e-4"<< endl;

// cin >> d;
cout << "Enter number of steps total: " << endl;
double stop = 0;
cin >> stop;
double po_per_cycle = 1e3; // again temporary
// cin >> po_per_cycle;
double mod_number = steps_per_cycle / po_per_cycle;
cout << "Enter data output filenname. Remember to add the extension: .txt!!!!"<< endl;
string out_filename;
cin >> out_filename;

ofstream ouch;
ouch.open(out_filename.c_str());
if (!ouch.good()) {
cerr << "Failed to open output file '" << out_filename << "'" << endl;
exit(1);
}
// column headers:
ouch << setw(15) << Q+"t in s"+Q << ","
<< setw(15) << Q+"theta in rad"+Q << ","
<< setw(15) << Q+"thetadot in rad/s"+Q << ","
<< setw(15) << Q+"deltag_now"+Q << ","
<< setw(15) << Q+"drive amp"+Q << ","
<< setw(15) << Q+"dissipation amp"+Q << ","
<< setw(15) << Q+"energy"+Q << ","
<< endl;

// print the initial state:
p(ouch, t, theta, thetadot, 0/*deltag*/, dr, r, energy);


int M = 0;

double deltag_now = 0;


//////////////////////////////////////////////////
// main loop
for (int N = 1; N <= stop; N++){


// zero crossing detection "previous" angle

double thetaprevious = theta;

//ODE solution
//Cromer solution:

double A = - ((4*pi*pi)/(period*period))* sin(theta) - r*thetadot*abs(thetadot) + dr;
thetadot = A * deltat + thetadot;
theta = theta + deltat * thetadot;

t = t + deltat; //step t for next iteration.
// end solution

/*
double A = - r*thetadot*abs(thetadot) +1 - ((4*pi*pi)/(period*period)) * sin(theta);

thetadot = A * deltat/2 + thetadot;

theta = theta + deltat * thetadot;

// recalculate acceleration at end of interval, using new theta:

A = - r*thetadot*abs(thetadot) +1 - ((4*pi*pi)/(period*period)) * sin(theta);

thetadot = thetadot + A*deltat/2;

t = t + deltat; // keep track of elapsed time
*/


//crossing detection continued

int M ;
double crossed = theta * thetaprevious;



if (crossed <= 0.0 ) //&& thetadot <= 0.0) {
{ M++;
}

/*
////////////////// //alternate method
float F = M + 0.5;

if (fmod (F, 15.5) == 0) {r=0;

} else {r=2;

}


////////////////
*/

energy = crossed;
int pulse = M; energy = pulse;


if ( pulse % 15 == 0) {deltag_now = 0;

} else { deltag_now = 1;

}


// above not used (delta_now == 0 missing from the "if" next line

if ( theta > -0.03 && theta < +0.03 && thetadot <0.0 ) { dr = -d;} else {dr = 0;}


pulse = energy;


// Print out something every mod_number steps, (or as close thereto as we can, if mod_number is not an integer).
// Also print out the very last result, (which will not be evenly spaced if total duration
// is not an integer multiple of the spacing between printouts).
if (N == stop || (mod(N + 0.5, mod_number)) < 1) {
p(ouch, t, theta, thetadot, deltag_now, dr, r, energy);
}
} // end main loop
}






//crossing detection continued

int M ;
double crossed = theta * thetaprevious;



if (crossed <= 0.0 ) //&& thetadot <= 0.0) {
{ M++;
}

float F = M + 0.5;

energy = crossed;
int pulse = M; energy = pulse;

if ( pulse % 10 == 0) {deltag_now = 0;

} else { deltag_now = 1;

}

if (deltag_now == 0 && theta > -0.03 && theta < +0.03 && thetadot <0.0 ) { dr = -d;} else {dr = 0;}

The problem w/ this code is the drive (dr is only one side of when theta is zero. To simulate the clock the drive must be "on" equally before and after when the pendulum is at bottom dead center, that is when the angle theta is zero. I'm too new w/ this type coding to solve this after several days work. I need it to submit to a horological science newsletter.

bgcm