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

1 Ansicht (letzte 30 Tage)
dcydhb dcydhb
dcydhb dcydhb am 11 Jun. 2019
Kommentiert: Jan am 12 Jun. 2019
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 Kommentare
Rik
Rik am 11 Jun. 2019
Bearbeitet: Rik am 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.
dcydhb dcydhb
dcydhb dcydhb am 12 Jun. 2019
thanks a lot!!!

Melden Sie sich an, um zu kommentieren.

Akzeptierte Antwort

Jan
Jan am 11 Jun. 2019
Bearbeitet: Jan am 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 Kommentare
dcydhb dcydhb
dcydhb dcydhb am 12 Jun. 2019
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
Jan am 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

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Kategorien

Mehr zu Call C++ from MATLAB finden Sie in Help Center und File Exchange

Produkte


Version

R2014a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by