gen_calc_model.pl 3.16 KB
Newer Older
jben's avatar
jben committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74
#!/usr/bin/perl

use strict;

sub maxima2cc
{
    my $chaine=$_[0];
    my $nomvariable=$_[1];
    my $indent=$_[2];
    my $provprefix=$_[3];
    my $old;
    my $out="";
    # number -> as double
    do
    {
        $old=$chaine;
        $chaine =~ s/([^A-Za-z_0-9\.]|^)([0-9]+)([^0-9\.]|$)/\1\2\.0\3/g;
    } while(!($old eq $chaine));

    # On recherche toutes les puissances pour les precalculer
    my %pui=();
    my $copy = $chaine;
    while($copy =~ s/([A-Za-z_0-9\.]+\^[0-9\.]+)//)
    {
        $pui{$1}++;
    }

    my $compteur=0;
    foreach my $expr (keys(%pui))
    {
        $expr=~m/(.*)\^(.*)/;
        $out="$out${indent}double $provprefix$compteur = std::pow($1,$2); // $expr\n";
        $expr =~ s/\^/\\\^/;
        $chaine =~ s/$expr/$provprefix$compteur/g;

        $compteur++;
    }


    # std::pow
    my $chaineside;
    $chaineside = qr/
            (
                (?>[A-Za-z_0-9\.]+)
            |
                \(
                    (?>[^()]*)
                    (
                        (?1)(?>[^()]*)
                    )*
                \)
            )
    /x;


    do
    {
        $old=$chaine;
        $chaine =~ s/(?<a>$chaineside)\^(?<b>$chaineside)/pow($+{'a'},$+{'b'})/g;
    } while(!($old eq $chaine));

    $chaine =~ s/pow/std::pow/g;

    # std::sqrt
    $chaine =~ s/sqrt/std::sqrt/g;

    # espacement
    $chaine =~ s/([\+\-\*\/])/ \1 /g;

    # fin
    $out= "$out\n$indent$nomvariable=$chaine\n\n";
    return $out;
}

75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90
sub get_result
{
    my $filename = $_[0];
    my $varname = $_[1];
    my $indent = $_[2];
    open(my $g2,"<$filename");
    my $inEa="";
    while(<$g2>)
    {
        chomp;
        $inEa="$inEa$_";
    }
    close($g2);
    my $r=maxima2cc($inEa,$varname,$indent,"d_prov_${varname}_");
    return $r;
}
jben's avatar
jben committed
91

92
# processing du calc_model.mc
jben's avatar
jben committed
93 94 95 96 97 98 99 100 101 102
system("maxima -b calc_model.mc");


# Lecture du resultat

# On ouvre la sortie
open(my $h,">calc_model.cc");
print $h "#include<cmath>\n\n";
print $h "void model::calc_model(double & d_V1, double & d_Ea, double & d_Ec, double d_V0, double d_vent, double d_d, double d_Ep, double d_Er) const\n";
print $h "{\n";
103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122
print $h "    double a,b,c,d,p,q,discri;\n\n";

# on calcul a,b,c,d puis p,q puis discri
print $h get_result("calc_model_a.out","a","    ");
print $h get_result("calc_model_b.out","b","    ");
print $h get_result("calc_model_c.out","c","    ");
print $h get_result("calc_model_d.out","d","    ");
print $h get_result("calc_model_p.out","p","    ");
print $h get_result("calc_model_q.out","q","    ");
print $h get_result("calc_model_discri.out","discri","    ");

# puis si discri > 0 on prend une solution, sinon l'autre
print $h "    if(discri>0)\n";
print $h "    {\n";
print $h get_result("calc_model_sol_pos.out","d_V1","        ");
print $h "    }\n";
print $h "    else\n";
print $h "    {\n";
print $h get_result("calc_model_sol_neg.out","d_V1","        ");
print $h "    }\n";
jben's avatar
jben committed
123 124 125 126 127 128 129 130


print $h "    if(d_V1<0)\n";
print $h "        d_V1=0;\n\n";

print $h "    if(std::isnan(d_V1))\n";
print $h "            d_V1=d_V0;\n\n";

131
# puis on calcule les energie
jben's avatar
jben committed
132

133 134
print $h get_result("calc_model_Ec.out","d_Ea","    ");
print $h get_result("calc_model_Ea.out","d_Ec","    ");
jben's avatar
jben committed
135 136 137

print $h "}\n";
close($h);