/* Included only by an isolated runtime copy. No production RHS mode, force
 * formula, state dimension, Jacobian coloring or property cache is changed. */
#include <float.h>

#if NEXPERIMENT_CONTACTS
static double experiment_penetration(int j,const double *y) {
    ExperimentContact c=experiment_contacts[j];
    if(c.kind==1) return c.boundary-y[c.v1+1];
    if(c.kind==2) return y[c.v1+1]-c.boundary;
    /* Preserve exactly the production gap expression and operation order. */
    return -(c.boundary+y[c.v2+1]-y[c.v1+1]);
}
static double experiment_velocity(int j,const double *y) {
    ExperimentContact c=experiment_contacts[j];
    if(c.kind==1) return -y[c.v1];
    if(c.kind==2) return y[c.v1];
    return y[c.v1]-y[c.v2];
}
static double experiment_value(int j,int force,const double *y) {
    double p=experiment_penetration(j,y);
    if(force==2) return experiment_velocity(j,y);
    if(!force) return p;
    ExperimentContact c=experiment_contacts[j];
    double fraction=c.pdis>0 ? -expm1(-fmax(p,0)/c.pdis) : 1;
    return c.stiffness*p+fraction*c.damping*experiment_velocity(j,y);
}
static double experiment_locate(NativeRun *r,int j,int force,double left,double right,
                                double sign,NativeDense dense,void *context) {
    double y[NSTATES];
    for(int k=0;k<60;k++) {
        double mid=left+.5*(right-left);
        if(mid<=left || mid>=right) break;
        r->experiment_dense++;r->experiment_roots++;
        if(!dense(context,mid,y)) return NAN;
        double g=experiment_value(j,force,y);
        if(!isfinite(g)) return NAN;
        if(sign<0 ? g>=0 : g<=0) right=mid; else left=mid;
    }
    return right;
}
static double experiment_bracket(NativeRun *r,int j,int force,double left,double right,
                                double a,double b,NativeDense dense,void *context,int *direction) {
    int slot=2*j+force;
    if(r->experiment_direction[slot] && r->experiment_last[slot]==left)
        a=r->experiment_direction[slot]*DBL_MIN;
    if(a==0 || (a>0 ? b>0 : b<0)) return INFINITY;
    *direction=a<0 ? 1 : -1;
    return experiment_locate(r,j,force,left,right,a,dense,context);
}
static double experiment_contact_candidate(NativeRun *r,int j,int force,double t,double next,
                                           const double *old,const double *trial,
                                           NativeDense dense,void *context,int *direction) {
    double a=experiment_value(j,force,old), b=experiment_value(j,force,trial);
    if(!isfinite(a) || !isfinite(b)) return NAN;
    /* Detect a gap excursion and return through the same boundary even when
     * the endpoint gaps have the same sign: split at the velocity reversal.
     * Like the existing event locator, this relies on resolved accepted steps;
     * arbitrarily many unresolved oscillations in one step are not certified. */
    if(!force) {
        double va=experiment_velocity(j,old),vb=experiment_velocity(j,trial);
        if((va<0 && vb>0) || (va>0 && vb<0)) {
            double turn=experiment_locate(r,j,2,t,next,va,dense,context),y[NSTATES];
            if(!isfinite(turn)) return NAN;
            r->experiment_dense++;
            if(!dense(context,turn,y)) return NAN;
            double g=experiment_value(j,force,y);
            if(g==0 && ((a<0 && b<0) || (a>0 && b>0))) return INFINITY;
            double found=experiment_bracket(r,j,force,t,turn,a,g,dense,context,direction);
            if(isfinite(found) || isnan(found)) return found;
            return experiment_bracket(r,j,force,turn,next,g,b,dense,context,direction);
        }
    }
    return experiment_bracket(r,j,force,t,next,a,b,dense,context,direction);
}
#endif

#if NEXPERIMENT_RELEASES
static int experiment_drives(NativeRun *r,double t,const double *y,double *drives) {
    r->nfev++;r->experiment_rhs++;
    return model_experiment_release_drives(t,y,drives);
}
static double experiment_release_root(NativeRun *r,int j,int lower,double left,double right,
                                      NativeDense dense,void *context) {
    double y[NSTATES],drives[NEXPERIMENT_RELEASES];
    for(int k=0;k<60;k++) {
        double mid=left+.5*(right-left);
        if(mid<=left || mid>=right) break;
        r->experiment_dense++;r->experiment_roots++;
        if(!dense(context,mid,y) || !experiment_drives(r,mid,y,drives)) return NAN;
        if(lower ? drives[j]>0 : drives[j]<0) right=mid; else left=mid;
    }
    return right;
}
#endif

static int experiment_candidates(NativeRun *r,double t,double next,const double *old,const double *trial,
                                 NativeDense dense,void *context,double *when,int *indices,int *kinds,int *count) {
    (void)r;(void)t;(void)next;(void)old;(void)trial;(void)dense;(void)context;
    (void)when;(void)indices;(void)kinds;(void)count;
#if NEXPERIMENT_CONTACTS
    for(int j=0;j<NEXPERIMENT_CONTACTS;j++) {
        r->experiment_checks++;
        for(int force=0;force<2;force++) {
            if(force && (experiment_contacts[j].signed_force==1 ||
                         experiment_penetration(j,old)<=0)) continue;
            int direction=0;
            double at=experiment_contact_candidate(r,j,force,t,next,old,trial,dense,context,&direction);
            if(isnan(at)) return 0;
            if(isfinite(at)) {
                if(force) {
                    double y[NSTATES];r->experiment_dense++;
                    if(!dense(context,at,y)) return 0;
                    if(experiment_penetration(j,y)<=0) continue;
                }
                int n=(*count)++;
                when[n]=at;indices[n]=j;kinds[n]=force?-3:-2;
                r->experiment_pending[2*j+force]=direction;
            }
        }
    }
#endif
#if NEXPERIMENT_RELEASES
    double before[NEXPERIMENT_RELEASES],after[NEXPERIMENT_RELEASES];
    int evaluated=0;
    for(int j=0;j<NEXPERIMENT_RELEASES;j++) {
        NativeStop s=model_stops[j];int v=s.velocity_index, x=v+1;
        if(fabs(old[v])>1e-12*fmax(fabs(old[v]),1)) continue;
        if(r->experiment_direction[2*NEXPERIMENT_CONTACTS+j] &&
           r->experiment_last[2*NEXPERIMENT_CONTACTS+j]==t) continue;
        for(int lower=0;lower<2;lower++) {
            double bound=lower?s.lower:s.upper,tol=1e-12*fmax(fabs(bound),1);
            if(lower ? old[x]>bound+tol : old[x]<bound-tol) continue;
            if(!evaluated) {
                if(!experiment_drives(r,t,old,before) || !experiment_drives(r,next,trial,after)) return 0;
                evaluated=1;
            }
            r->experiment_checks++;
            if(lower ? !(before[j]<=0 && after[j]>0) : !(before[j]>=0 && after[j]<0)) continue;
            double at=experiment_release_root(r,j,lower,t,next,dense,context);
            if(!isfinite(at)) return 0;
            int n=(*count)++;when[n]=at;indices[n]=j;kinds[n]=lower?-4:-5;
        }
    }
#endif
    return 1;
}

static void experiment_commit(NativeRun *r,double t,double *y,const double *when,const int *indices,const int *kinds,int count) {
    (void)y;
    for(int i=0;i<count;i++) {
        if(kinds[i]>-2 || when[i]!=t) continue;
        int j=indices[i];
        (void)j;
#if NEXPERIMENT_CONTACTS
        if(kinds[i]==-2 || kinds[i]==-3) {
            int force=kinds[i]==-3,slot=2*j+force;
            r->experiment_last[slot]=t;
            r->experiment_direction[slot]=r->experiment_pending[slot];
            if(force) r->experiment_clipping++;
            else if(experiment_contacts[j].kind==3) r->experiment_lstp++;
            else r->experiment_mass++;
            fprintf(stderr,"{\"phase\":\"mechanical-event\",\"kind\":\"%s\",\"index\":%d,\"time\":%.17g,\"direction\":%d,\"penetration\":%.17g,\"relativeVelocity\":%.17g}\n",
                force?"force-clip":experiment_contacts[j].kind==3?"lstp-contact":"mass-contact",j,t,
                r->experiment_direction[slot],experiment_penetration(j,y),experiment_velocity(j,y));
        }
#endif
#if NEXPERIMENT_RELEASES
        if(kinds[i]==-4 || kinds[i]==-5) {
            NativeStop s=model_stops[j];int slot=2*NEXPERIMENT_CONTACTS+j;
            y[s.velocity_index]=0;y[s.velocity_index+1]=kinds[i]==-4?s.lower:s.upper;
            r->experiment_last[slot]=t;r->experiment_direction[slot]=1;r->experiment_release++;
            fprintf(stderr,"{\"phase\":\"mechanical-event\",\"kind\":\"mass-release\",\"index\":%d,\"time\":%.17g,\"side\":\"%s\"}\n",
                j,t,kinds[i]==-4?"lower":"upper");
        }
#endif
    }
    (void)r;
}