How to use matlab to call c + function file and draw?
1 view (last 30 days)
Show older comments
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](https://www.mathworks.com/matlabcentral/answers/uploaded_files/223892/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
Rik
on 11 Jun 2019
Edited: Rik
on 11 Jun 2019
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.
Accepted Answer
Jan
on 11 Jun 2019
Edited: Jan
on 11 Jun 2019
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
Jan
on 12 Jun 2019
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
More Answers (0)
See Also
Categories
Find more on Call C++ from MATLAB in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!