#include "stdafx.h"

#include "DllFunc.h"

#include "math.h"

TireForce_API void __cdecl tire_force

           (double time, double itime, double conprop[], double tireprop[], double massprop[],

           double upar[], int npar, int jflag, int iflag, double tfsae[3], double rfsae[3], double fcoef[2])

{

           using namespace rd_syscall;

           // Parameter Information

           //   time     : Simulation time of RD/Solver. (Input)

           //   itime    : Past time for saved values. (Input)

           //   conprop  : Contact property for computung Tire force. (Input)

           //   tireprop : Geometry and Material property of TIRE. (Input)

           //   massprop : Mass property of TIRE. (Input)

           //   upar     : Parameters defined by user. (Input)

           //   npar     : Number of user parameters. (Input)

           //   jflag    : When RD/Solver evaluates a Jacobian, the flag is true. (Input)

           //   iflag    : When RD/Solver initializes arrays, the flag is true. (Input)

           //   tfsae    : Returned froce vector acting at the contact point in the SAE coordinate system. (Output, Size : 3)

           //   rfsae    : Returned torque vector acting at the contact point in the SAE coordinate system. (Output, Size : 3)

           //   fcoef    : Returned friction coefficients. (Output, Size : 2)

           // User Statement

           double value;

           double zkForce, zcForce, new_zcForce;

           double ver_stiff, pen, mass, damping, dot_pen, max_pen;

           int errflg;

           //step function

           pen = conprop[3];

           max_pen = 0.01* tireprop[1];

           rd_step(pen, 0, 0, max_pen, 1, 1, &value, &errflg);

           //calculation for normal force

           ver_stiff = tireprop[2];

           zkForce = -ver_stiff*pen;

           mass = massprop[0];

           damping = tireprop[7];

           dot_pen = conprop[4];

           new_zcForce = -2.0*sqrt(mass*abs(ver_stiff))*damping*dot_pen;

           zcForce = value*new_zcForce;

           //z axis normal force

           tfsae[2] = min(0.0, zkForce+zcForce);

}