MATLAB Answers

0

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

Asked by dcydhb dcydhb on 11 Jun 2019
Latest activity Commented on by Jan
on 12 Jun 2019
Accepted Answer by Jan
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?
cd-kc.png
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;
}

  4 Comments

Show 1 older comment
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.
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.

Sign in to comment.

Products


Release

R2014a

1 Answer

Answer by Jan
on 11 Jun 2019
Edited by Jan
on 11 Jun 2019
 Accepted Answer

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 .

  3 Comments

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.
so in fact i need your answer to call the c++ in the matlab or just use the c++ to call rather than rewrite.
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

Sign in to comment.