#include "stdafx.h"

#include "DllFunc.h"

#include "math.h"

ModalForceExt_API void __cdecl modal_force_ext

           (int id, double time, double upar[], int npar, int ifbody, int nodarr[], int nonde, int jflag, int iflag, double result[])

{

           using namespace rd_syscall;

           // Parameter Information

           //   id     : Modal force sequential identification. (Input)

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

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

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

           //   ifbody : RFLEX Body sequential ID. (Input)

           //   nodarr : Node ID array of input node set. (Input)

           //   nonde  : Number of node of node set. (Input)

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

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

           //   result : Returned nodal force vector. Acting point of the nodal force is each center of each node.

           //            Reference frame of each force vector must be the ground inertia marker. (Output, Size: nonde * 6)

           // User Statement

           int i, j, errflg;

           int nodeseq;

           double Area, Cd, rho;

           double vn[3], vn_mag, ucf, scale;

           WCHAR NodeName[256];

           char  strNodeName[256];

           Cd = upar[0];

           rho = upar[1];

           scale = upar[2];

           rd_ucf(&ucf);

           for(i=0;i<nonde;i++){

                     if(nodarr[i] == 50004 || nodarr[i] == 50007) Area = 100;

                     else if (nodarr[i] >=50010 && nodarr[i] <= 50075 ) Area = 200;

                     else if(nodarr[i] > 50075) Area =400;

                     // Set velocity vector of node

                     get_rflex_nodeseqid(ifbody,nodarr[i],&nodeseq,&errflg);

                     get_rflex_nodetvel(ifbody,nodeseq,NULL,0,vn,&errflg);

                     vn_mag = sqrt(pow(vn[0],2)+pow(vn[1],2)+pow(vn[2],2));

                     for(j=0;j<3;j++){

                                if(vn_mag == 0.0)

                                           result[i*6+j] = 0.0;

                                else

                                           result[i*6+j] = -1.0/2.0*Cd*rho*vn_mag*vn_mag*Area

                                           *(vn[j]/vn_mag) // Set negative direction of the nodal velocity

                                           /ucf            // Making force unit

                                           *scale;         //

                                result[i*6+j+3] = 0.0;

                     }

           }

}