How can I get coder to generate c code for a script with the solve function

I am trying to get c code from a Matlab algorithm using Matlab coder but the code uses the solve function which is not supported in coder. I tried solving for the variable that I am trying to obtain but I do not get the same answer. I trying to solve for current in the equation below:

                syms current
                eqn = log((pi()*exp(1)*(radius(l).^2)*CD)/deltaE_max) == log(current);
                soln2 = solve(eqn, current);
                I_cath(l) = soln2;

Can anyone tell me how to solve for current without using the solve function?

[Markus Hahnenkamm]

I am currently experiencing the same issue. I would like to integrate the solve in a c code. Is there a workaround?

[Walter Roberson]

Markus Hahnenkamm:

There is no way to include solve() in generated code. You need to attempt to find a closed form solution in MATLAB and then code that solution to get a numeric result. You might need to use fzero() or fsolve(). You will probably find it useful to look at the 'file' option of matlabFunction()

Answer 1

[Walter Roberson]

Your code is

eqn = log((pi()*exp(1)*(radius(l).^2)*CD)/deltaE_max) == log(current);

Current is only found on the right hand side, in a fairly simple usage. To solve for it, take exp() of both sides:

exp(log((pi()*exp(1)*(radius(l).^2)*CD)/deltaE_max)) == exp(log(current))

which gives

(pi*CD*exp(1)*radius(l)^2)/deltaE_max == current

which gives a direct calculation for current based upon the other values. solve() is not needed.

Answer 2

[Cordell Smith]

Sorry I copied the equation incorrectly here is the actually equation:

          for l = 1:radius_el
            syms current
            eqn = (4*pi()*A*B*deltaE_max) / current + log((pi()*exp(1)*(radius(l).^2)*CD)/deltaE_max) == log(current);
            soln = solve(eqn, current);

current is actually on both sides of the equation so it was not such a simple solution. I used a work around by setting both sides of the equation equal to each other and got the approximate solutions by finding where they intersected using the code below:

for l = 1:radius_el

            current=[.00001:0.00000001:.01];
         y2=log(current);
          y3=(4*pi()*A*B*deltaE_max) ./ current + log((pi()*exp(1)*(radius(l).^2)*CD)/deltaE_max);
           dif=abs(y2-y3);

             soln=find(dif==min(dif(:)))+1000;

[Walter Roberson]

syms A B deltaE_max radius(l) CD 
eqn = (4*sym('pi')*A*B*deltaE_max) / current + log((sym('pi')*exp(sym(1))*(radius(l).^2)*CD)/deltaE_max) == log(current);
soln = solve(eqn, current)

This gives

(4*A*B*deltaE_max*pi)/wrightOmega(- log(1/(4*A*B*deltaE_max*pi)) - log((pi*CD*exp(1)*radius(l)^2)/deltaE_max))

Most of this can have code generated. The part that cannot is the WrightOmega. It is probably easier to find a replacement code if you rewrite in terms of LambertW, which happens to give

4*A*B*deltaE_max*pi / lambertw(4*A*B*deltaE_max^2*exp(-1)/(CD*radius(l)^2))

In that regard, I find https://github.com/DarkoVeberic/LambertW . Or you might be satisfied with using http://lasp.colorado.edu/cism/CISM_DX/code/CISM_DX-0.50/required_packages/octave-forge/main/specfun/lambertw.m

[Walter Roberson]

In the special case that A, B, deltaE_max, and CD are fixed ahead of time, and radius does not vary much, you could consider doing a taylor expansion. This would involve a number of uses of LambertW(4*A*B*deltaE_max^2*exp(-1)/(nominal_R^2*CD)) but under the special situation I mention that could be calculated ahead of time.