How to use matlab to call c + function file and draw？

The following c++ function file is a formula that contains some piecewise functions. Let's first assume that Re=2320, so how to uses matlab to call this function file to make a curve of cd versus Kc.

The range of Kc is 0 to 100.

i just transform many codes and use the matlab codes to plot and get the figure,but can we just call c++ codes and plot in the matlab?

the c++ function file is as this

#include"force.h"
#include"lexer.h"
#include"fdm.h"
#include"ghostcell.h"

double force::morison_cd(lexer* p, fdm *a, ghostcell *pgc, double Re, double Kc)
{
    cd=0.0;
    double cd6,cd8,cd10,cd15,cd20,cd40,cd60,cd100;

    if(Kc<=8.0)
    {
    cd6 = - 5.373e-51*pow(Re,9.0)
          + 2.486e-44*pow(Re,8.0)
          - 4.879e-38*pow(Re,7.0)
          + 5.300e-32*pow(Re,6.0)
          - 3.941e-26*pow(Re,5.0)
          + 1.437e-20*pow(Re,4.0)
          - 3.677e-15*pow(Re,3.0)
          + 5.605e-10*pow(Re,2.0)
          - 4.507e-06*Re
          + 1.913;
    cd6 = MIN(cd6,1.5);
    }

    if(Kc<=10.0 && Kc>6.0)
    {
    cd8 = - 3.111e-51*pow(Re,9.0)
          + 1.461e-44*pow(Re,8.0)
          - 2.947e-38*pow(Re,7.0)
          + 3.340e-32*pow(Re,6.0)
          - 2.336e-26*pow(Re,5.0)
          + 1.040e-20*pow(Re,4.0)
          - 2.927e-15*pow(Re,3.0)
          + 4.978e-10*pow(Re,2.0)
          - 4.491e-05*Re
          + 2.203;
    cd8 = MIN(cd8,1.765);
    }

    if(Kc<=15.0 && Kc>8.0)
    {
    cd10 = - 6.845e-51*pow(Re,9.0)
           + 3.125e-44*pow(Re,8.0)
           - 6.067e-38*pow(Re,7.0)
           + 6.540e-32*pow(Re,6.0)
           - 4.295e-26*pow(Re,5.0)
           + 1.775e-20*pow(Re,4.0)
           - 4.605e-15*pow(Re,3.0)
           + 7.212e-10*pow(Re,2.0)
           - 6.037e-05*Re
           + 2.69;
    cd10 = MIN(cd10,2.0);
    }

    if(Kc<=20.0 && Kc>10.0)
    {
    cd15 = - 2.991e-51*pow(Re,9.0)
           + 1.506e-44*pow(Re,8.0)
           - 3.222e-38*pow(Re,7.0)
           + 3.824e-32*pow(Re,6.0)
           - 2.755e-26*pow(Re,5.0)
           + 1.243e-20*pow(Re,4.0)
           - 3.505e-15*pow(Re,3.0)
           + 5.985e-10*pow(Re,2.0)
           - 5.585e-05*Re
           + 2.804;
    cd15 = MIN(cd15,2.0);
    }

    if(Kc<=40.0 && Kc>15.0)
    {
    cd20 = - 1.184e-51*pow(Re,9.0)
           + 6.604e-45*pow(Re,8.0)
           - 1.561e-38*pow(Re,7.0)
           + 2.044e-32*pow(Re,6.0)
           - 1.625e-26*pow(Re,5.0)
           + 8.103e-21*pow(Re,4.0)
           - 2.516e-15*pow(Re,3.0)
           + 4.660e-10*pow(Re,2.0)
           - 4.607e-05*Re
           + 2.481;
    cd20 = MIN(cd20,1.96);
    }

    if(Kc<=60.0 && Kc>20.0)
    cd40 = 2.480e-46*pow(Re,8.0)
         - 1.277e-39*pow(Re,7.0)
         + 2.762e-33*pow(Re,6.0)
         - 3.263e-27*pow(Re,5.0)
         + 2.290e-21*pow(Re,4.0)
         - 9.706e-16*pow(Re,3.0)
         + 2.406e-10*pow(Re,2.0)
         - 3.133e-05*Re
         + 2.185;

    if(Kc<=100.0 && Kc>40.0)
    {
    cd60 = 3.705e-46*pow(Re,8.0)
         - 1.644e-39*pow(Re,7.0)
         + 3.065e-33*pow(Re,6.0)
         - 3.133e-27*pow(Re,5.0)
         + 1.933e-21*pow(Re,4.0)
         - 7.503e-16*pow(Re,3.0)
         + 1.841e-10*pow(Re,2.0)
         - 2.644e-05*Re
         + 2.187;
    cd60 = MIN(cd60,1.8);
    }

    if(Kc>60.0)
    {
    cd100 = - 1.430e-51*pow(Re,9.0)
            + 6.290e-45*pow(Re,8.0)
            - 1.171e-38*pow(Re,7.0)
            + 1.204e-32*pow(Re,6.0)
            - 7.533e-27*pow(Re,5.0)
            + 3.004e-21*pow(Re,4.0)
            - 7.962e-16*pow(Re,3.0)
            + 1.507e-10*pow(Re,2.0)
            - 2.096e-05*Re
            + 2.097;
    cd100 = MIN(cd100,1.4);
    }




    if(Kc<=6.0)
    cd = cd6;

    if(Kc<=8.0 && Kc>6.0)
    cd = ((8.0-Kc)/2.0)*cd6 + ((Kc-6.0)/2.0)*cd8;

    if(Kc<=10.0 && Kc>8.0)
    cd = ((10.0-Kc)/2.0)*cd8 + ((Kc-8.0)/2.0)*cd10;

    if(Kc<=15.0 && Kc>10.0)
    cd = ((15.0-Kc)/5.0)*cd10 + ((Kc-10.0)/5.0)*cd15;

    if(Kc<=20.0 && Kc>15.0)
    cd = ((20.0-Kc)/5.0)*cd15 + ((Kc-15.0)/5.0)*cd20;

    if(Kc<=40.0 && Kc>20.0)
    cd = ((40.0-Kc)/20.0)*cd20 + ((Kc-20.0)/20.0)*cd40;

    if(Kc<=60.0 && Kc>40.0)
    cd = ((60.0-Kc)/20.0)*cd40 + ((Kc-40.0)/20.0)*cd60;

    if(Kc<=100.0 && Kc>60.0)
    cd = ((100.0-Kc)/40.0)*cd60 + ((Kc-60.0)/40.0)*cd100;

    if(Kc>100.0)
    cd = cd100;


    return cd;

}

[Rik]

You could compile the code to a mex function. But I would suggest rewriting your function in Matlab, so you can take advantage of the array processing capabilities. That way you can avoid using a loop.

[dcydhb dcydhb]

do you mean that if i use call the c++ code i must use a loop,why?

in fact when rewrite in the matlab,there are many item that you have to change.

[Rik]

I have very limited experience with C++, so if my following assumption is false you can ignore this. It looks like you can only pass scalars as an input to your function. So if you want to use this function to plot a line, you will have to call the function in some sort of a loop to retrieve the y-value for each x-value separately.

If you re-work this to use logical indexing you can take advantage of array processing, which should result in a considerable speed boost compared to a loop, even when using a mex function.

[dcydhb dcydhb]

thanks a lot!!!

Answer 1

[Jan]

What about converting the code to Matlab? This is more or less trivial:

function cdValue = yourFcn(Kc, Re)
cdValue=0.0;  % Do not shadow Matlab's CD function
if Kc<=8.0
    cd6 = - 5.373e-51*pow(Re,9.0) ...
          + 2.486e-44*pow(Re,8.0) ...
          - 4.879e-38*pow(Re,7.0) ...
          + 5.300e-32*pow(Re,6.0) ...
          - 3.941e-26*pow(Re,5.0) ...
          + 1.437e-20*pow(Re,4.0) ...
          - 3.677e-15*pow(Re,3.0) ...
          + 5.605e-10*power(Re,2.0) ...
          - 4.507e-06*Re ...
          + 1.913;
    cd6 = min(cd6,1.5);
elseif Kc<=10.0 && Kc>6.0
    cd8 = - 3.111e-51*pow(Re,9.0) ...
          + 1.461e-44*pow(Re,8.0) ...
          - 2.947e-38*pow(Re,7.0) ...
          + 3.340e-32*pow(Re,6.0) ...
          - 2.336e-26*pow(Re,5.0) ...
          + 1.040e-20*pow(Re,4.0) ...
          - 2.927e-15*pow(Re,3.0) ...
          + 4.978e-10*pow(Re,2.0) ...
          - 4.491e-05*Re ...
          + 2.203; ...
    cd8 = min(cd8,1.765);
elseif ... and so on

if Kc <= 6.0
    cdValue = cd6;
elseif(Kc<=8.0 && Kc>6.0)
    cdValue = ((8.0-Kc)/2.0)*cd6 + ((Kc-6.0)/2.0)*cd8;
elseif(Kc<=10.0 && Kc>8.0)
    cdValue = ((10.0-Kc)/2.0)*cd8 + ((Kc-8.0)/2.0)*cd10;
elseif(Kc<=15.0 && Kc>10.0)
    cdValue = ((15.0-Kc)/5.0)*cd10 + ((Kc-10.0)/5.0)*cd15;
elseif(Kc<=20.0 && Kc>15.0)
    cdValue = ((20.0-Kc)/5.0)*cd15 + ((Kc-15.0)/5.0)*cd20;
elseif(Kc<=40.0 && Kc>20.0)
    cdValue = ((40.0-Kc)/20.0)*cd20 + ((Kc-20.0)/20.0)*cd40;
elseif(Kc<=60.0 && Kc>40.0)
    cdValue = ((60.0-Kc)/20.0)*cd40 + ((Kc-40.0)/20.0)*cd60;
elseif(Kc<=100.0 && Kc>60.0)
    cdValue = ((100.0-Kc)/40.0)*cd60 + ((Kc-60.0)/40.0)*cd100;
elseif(Kc>100.0)
    cdValue = cd100;
end

This took me 2 minutes now.

It would be smarter to vectorize the loop and to use logical indices instead of a pile of if elseif blocks. But if this function is not called millions of times, this version is fine also. Simply call it in a loop:

function main
        
KcVector =         
end

Sorry, I'd really wanted to write this loop for you also to show, how easy this is. Unfortunately the interface of this forum impedes typing code such massively that I'm really annoyed. I give up here after struggeling too hard with editing some lines of code. See https://www.mathworks.com/matlabcentral/answers/438664-strange-behavior-of-the-editor-in-the-forum .

[dcydhb dcydhb]

it is because the definition od the piecewise is Too chaotic and it has a lot of overlapping intervals so may be if you use the matlab to rewrite the piecewise function,you should rewrite and redefine the intervals.

so i wonder if there is the way just write a few lines in c++ or matlab to call it rather than rewrite the function files.

[dcydhb dcydhb]

so in fact i need your answer to call the c++ in the matlab or just use the c++ to call rather than rewrite.

[Jan]

Sorry, what? I've converted the posted function for you already. All you have to do is to write a for loop, which calls this function with the given inputs and store the output in a vector. I try it again hoping that the forum's interface let me do this:

function main
Re = 2320;
KcVec = 0:100;
cdVec = zeros(size(KcVec));
for k = 1:numel(KcVec)
    cdVec(k) = yourFcn(KcVec(k), Re);
end
end

This should be working. The only problem I had was the interface of the Matlab Answers forum, which applies too many "smart" guesses, re-indents my code repeatedly, inserts "end" statements, when it likes to do so, scrolls the section for editing out of the window, etc.

If you really want to run this through the mex interface, I'd prefer a C-function:

!!!UNTESTED CODE!!!

#include "mex.h"
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
double *Kc, Re, *cd, *KcEnd;

Kc    = mxGetPr(prhs[0]);
KcEnd = Kc + mxGetNumberOfElements(prhs[0]);
Re    = mxGetScalar(prhs[1]);

plhs[0] = mxCreateDoubleMatrix(1, mxGetNumberOfElements(prhs[0]), mxREAL);
cd      = mxGetPr(plhs[0]);

while (Kc < KcEnd) {
   *cd++ = morison_cd(Re, *Kc++);
}
return;
}

double morison_cd(double Re, double Kc)
{
    double cd=0.0;
    ... and your function

Compile this with:

    mex -R2017b yourFcn.c