%
%  arbit_force_function.m  ver 2.0  (batch-recursive implementation)
%  by Tom Irvine (modified)
%
function [acceleration,velocity,displacement,trans_force] = ...
    arbit_force_function(t,f,dt,mass,fn,damp,n_force_mass_unit)

    fig_num = 1;
    n = length(t);

    %--------------------------------------------------------------
    % Units handling
    %--------------------------------------------------------------
    if (n_force_mass_unit == 1)
        % input mass in lbm → convert to slug
        mass = mass/386;
    end

    fprintf('mass=%g\n',mass);

    %--------------------------------------------------------------
    % Initialize coefficients
    %--------------------------------------------------------------
    disp(' ');
    disp(' Initialize coefficients');

    omegan = 2.*pi*fn;
    omegad = omegan*sqrt(1 - damp^2);

    fprintf(' omegan=%g   omegad=%g \n',omegan,omegad);
   

    cosd = cos(omegad*dt);
    sind = sin(omegad*dt);

    fprintf(' cosd=%g   sind=%g \n',cosd,sind);
   

    domegadt = damp*omegan*dt;

    %--------------------------------------------------------------
    % Smallwood ramp-invariant coefficients (same as original)
    %--------------------------------------------------------------
    eee1  = exp(-domegadt);
    eee2  = exp(-2.*domegadt);

    ecosd = eee1*cosd;
    esind = eee1*sind;

    a1 =  2.*ecosd;
    a2 = -eee2;

    % sprintf(' a1=%g  a2=%g ',a1(1),a2(1));

    omeganT = omegan*dt;

    phi = 2*(damp)^2 - 1;
    DD1 = (omegan/omegad)*phi;
    DD2 = 2*DD1;

    df1 =  2*damp*(ecosd - 1) + DD1*esind + omeganT;
    df2 = -2*omeganT*ecosd    + 2*damp*(1 - eee2) - DD2*esind;
    df3 = (2*damp + omeganT)*eee2 + (DD1*esind - 2*damp*ecosd);

    VV1 = -(damp*omegan/omegad);

    vf1 = (-ecosd + VV1*esind) + 1;
    vf2 =  eee2 - 2*VV1*esind - 1;
    vf3 =  ecosd + VV1*esind  - eee2;

    MD  = (mass*omegan^3*dt);
    df1 = df1/MD;
    df2 = df2/MD;
    df3 = df3/MD;

    VD  = (mass*omegan^2*dt);
    vf1 = vf1/VD;
    vf2 = vf2/VD;
    vf3 = vf3/VD;

    af1 = esind/(mass*omegad*dt);
    af2 = -2*af1;
    af3 = af1;

    %--------------------------------------------------------------
    % Residual handling (unchanged – currently disabled)
    %--------------------------------------------------------------
    ire = 2;   % 1=yes, 2=no

    if (ire == 1)
        % Add trailing zeros for residual response
        disp(' ');
        disp(' Add trailing zeros for residual response');
        tmax  = (tmx - tmi) + 2./fn;
        limit = round(tmax/dt);
        n     = limit;
        % If you enable this, you must also pad f, t accordingly.
    end

    %--------------------------------------------------------------
    % Batch-recursive SDOF engine (displacement, velocity, acceleration)
    %--------------------------------------------------------------
    % Shared denominator for all three responses:
    %  y[n] = a1*y[n-1] + a2*y[n-2] + b1*x[n] + b2*x[n-1] + b3*x[n-2]
    %
    % Here, we batch three systems:
    %  Row 1: displacement  (df1, df2, df3)
    %  Row 2: velocity      (vf1, vf2, vf3)
    %  Row 3: acceleration  (af1, af2, af3)
    %
    B = [ df1, df2, df3;
          vf1, vf2, vf3;
          af1, af2, af3 ];

    f = f(:);              % ensure column
    N = length(f);

    % Pre-shifted force histories x(n), x(n-1), x(n-2)
    x0 = f;                            % x[n]
    x1 = [0; f(1:end-1)];              % x[n-1]
    x2 = [0; 0; f(1:end-2)];           % x[n-2]

    % Initialize output matrix: 3 rows (d,v,a), N columns (time steps)
    Y  = zeros(3,N);

    % Previous outputs for recursion
    y_prev1 = zeros(3,1);  % y[n-1]
    y_prev2 = zeros(3,1);  % y[n-2]

    % Main batch recursion loop
    for k = 1:N
        xin   = [x0(k); x1(k); x2(k)];      % [x(n); x(n-1); x(n-2)]
        y_new = a1*y_prev1 + a2*y_prev2 + B*xin;
        Y(:,k) = y_new;

        % Shift states
        y_prev2 = y_prev1;
        y_prev1 = y_new;
    end

    % Extract displacement, velocity, acceleration responses
    d_resp = Y(1,:).';  % column
    v_resp = Y(2,:).';
    a_resp = Y(3,:).';

    Q = 1/(2*damp);

    %--------------------------------------------------------------
    % Unit scaling for displacement (same as original)
    %--------------------------------------------------------------
    if (n_force_mass_unit == 2)
        % SI → displacement in mm if original convention
        d_resp = d_resp*1000;
    end

    %--------------------------------------------------------------
    % Package time histories
    %--------------------------------------------------------------
    displacement(:,1) = t(:);
    displacement(:,2) = d_resp;

    velocity(:,1)     = t(:);
    velocity(:,2)     = v_resp;

    acceleration(:,1) = t(:);
    acceleration(:,2) = a_resp;  % (will be converted to g below)

    %--------------------------------------------------------------
    % Transmitted force
    %--------------------------------------------------------------
    trans_force(:,1) = t(:);
    trans_force(:,2) = f - mass*acceleration(:,2);   % F_trans = F - m*a

    %--------------------------------------------------------------
    % Convert acceleration to g-units for plotting / output
    %--------------------------------------------------------------
    if (n_force_mass_unit == 1)
        a_resp_g = a_resp/386;
    else
        a_resp_g = a_resp/9.81;
    end

    acceleration(:,2) = a_resp_g;  % overwrite with g

    %--------------------------------------------------------------
    % Simple accel plot
    %--------------------------------------------------------------
    figure(fig_num); fig_num = fig_num + 1;
    plot(t,a_resp_g);
    xlabel('Time (sec)');
    ylabel('Accel (G)');
    if (n_force_mass_unit == 1)
        ttl = sprintf('SDOF Acceleration Response: mass=%g lbm  fn=%g Hz Q=%g',mass*386,fn,Q);
    else
        ttl = sprintf('SDOF Acceleration Response: mass=%g kg  fn=%g Hz Q=%g',mass,fn,Q);
    end
    title(ttl);
    grid on;
    zoom on;

    %--------------------------------------------------------------
    % Titles / labels for your multi-plot function
    %--------------------------------------------------------------
    tt = t(:);
    d  = d_resp;
    v  = v_resp;
    a  = a_resp_g;

    ay = 'Accel (G)';
    if (n_force_mass_unit == 1)
        aat = sprintf('SDOF Acceleration Response: mass=%g lbm  fn=%g Hz Q=%g',mass*386,fn,Q);
    else
        aat = sprintf('SDOF Acceleration Response: mass=%g kg  fn=%g Hz Q=%g',mass,fn,Q);
    end

    if (n_force_mass_unit == 1)
        vy  = 'Vel (inch/sec)';
        vvt = sprintf('SDOF Velocity Response: mass=%g lbm  fn=%g Hz Q=%g',mass*386,fn,Q);
    else
        vy  = 'Vel (m/sec)';
        vvt = sprintf('SDOF Velocity Response: mass=%g kg  fn=%g Hz Q=%g',mass,fn,Q);
    end

    if (n_force_mass_unit == 1)
        dy  = 'Disp (inch)';
        ddt = sprintf(' SDOF Displacement Response: mass=%g lbm fn=%g Hz Q=%g ',mass*386,fn,Q);
    else
        dy  = 'Disp (mm)';
        ddt = sprintf(' SDOF Displacement Response: mass=%g kg fn=%g Hz Q=%g ',mass,fn,Q);
    end

    [fig_num,~] = plot_avd_time_histories_subplots_altp_titles( ...
                             fig_num,tt,a,v,d,ay,vy,dy,aat,vvt,ddt);

    %--------------------------------------------------------------
    % Applied vs transmitted force subplot
    %--------------------------------------------------------------
    data1   = [t(:), f(:)];
    data2   = trans_force;

    xlabel2 = 'Time (sec)';

    if (n_force_mass_unit == 1)
        ylabel1 = 'Force (lbf)';
    else
        ylabel1 = 'Force (N)';
    end
    ylabel2 = ylabel1;

    t_string1 = 'Applied Force';
    if (n_force_mass_unit == 1)
        t_string2 = sprintf('Transmitted Force: mass=%g lbm  fn=%g Hz Q=%g',mass*386,fn,Q);
    else
        t_string2 = sprintf('Transmitted Force: mass=%g kg  fn=%g Hz Q=%g',mass,fn,Q);
    end

    fig_num = subplots_two_linlin_two_titles_scale_same( ...
                    fig_num,xlabel2,ylabel1,ylabel2, ...
                    data1,data2,t_string1,t_string2);

end
