load MuscleParaShortInvert.mat;
muscle = muscleInvert;
%muscle data structure including name, fmax, lceopt, lslack, width, order, dof, dofname, paranum, para
% length of femur 0.396m from SIMM
% length of tibia 0.43m from SIMM
Rmuscle = muscle;
for i=1:31
    Rmuscle(i).lceopt = muscle(i).lceopt * length(1,2)/0.396; %scaled by femur
    Rmuscle(i).lslack = muscle(i).lslack * length(1,2)/0.396; %scaled by femur
    [mrow, ncolumn] = size(muscle(i).para);
    paracolumn = ncolumn;
    Rmuscle(i).para(:,paracolumn) = muscle(i).para(:,paracolumn)*length(1,2)/0.396; %scaled by femur
end
for i=32:43
    Rmuscle(i).lceopt = muscle(i).lceopt * length(1,4)/0.43; %scaled by tibia
    Rmuscle(i).lslack = muscle(i).lslack * length(1,4)/0.43; %scaled by tibia
    [mrow, ncolumn] = size(muscle(i).para);
    paracolumn = ncolumn;
    Rmuscle(i).para(:,paracolumn) = muscle(i).para(:,paracolumn)*length(1,4)/0.43; %scaled by tibia
end
Lmuscle = muscle;
for i=1:31
    Lmuscle(i).lceopt = muscle(i).lceopt * length(1,2)/0.396; %scaled by femur
    Lmuscle(i).lslack = muscle(i).lslack * length(1,2)/0.396; %scaled by femur
    [mrow, ncolumn] = size(muscle(i).para);
    paracolumn = ncolumn;
    Lmuscle(i).para(:,paracolumn) = muscle(i).para(:,paracolumn)*length(1,2)/0.396; %scaled by femur
end
for i=32:43
    Lmuscle(i).lceopt = muscle(i).lceopt * length(1,4)/0.43; %scaled by tibia
    Lmuscle(i).lslack = muscle(i).lslack * length(1,4)/0.43; %scaled by tibia
    [mrow, ncolumn] = size(muscle(i).para);
    paracolumn = ncolumn;
    Lmuscle(i).para(:,paracolumn) = muscle(i).para(:,paracolumn)*length(1,4)/0.43; %scaled by tibia
end






fprintf(fid, 'muscle default \n');
fprintf(fid, '       a 0.25 \n');
fprintf(fid, '       vmrel 10 \n');
fprintf(fid, '       umax 0.04 \n');
fprintf(fid, '       fecmax 1.5 \n');
fprintf(fid, '       krel 1.0 \n');
fprintf(fid, '       slopfac 2.0 \n');
fprintf(fid, '       PEEslack 1.0 \n');
fprintf(fid, '       time_constants 0.050 0.060 \n');
fprintf(fid, '       end \n\n');

j =[1:10 12:16 20:24 28:35 38:40];
for k = 1:31
    i=j(k);
    fprintf(fid, 'muscle R%s \n', Rmuscle(i).name);
    fprintf(fid, '       fmax %6.1f \n', Rmuscle(i).fmax);
    fprintf(fid, '       lceopt %6.4f \n', Rmuscle(i).lceopt);
    fprintf(fid, '       lslack %6.4f \n', Rmuscle(i).lslack);
    fprintf(fid, '       width %5.3f \n', Rmuscle(i).width);
    if Rmuscle(i).peeslack>1.0
        fprintf(fid, '          PEEslack %5.3f \n', Rmuscle(i).peeslack);
    end
    fprintf(fid, '       geometry polynomial \n');
    fprintf(fid, '          dofs \n');
    fprintf(fid, '                R%s \n', Rmuscle(i).nest.dof1);
    if Rmuscle(i).nest.dof2 ~= []
        fprintf(fid, '                R%s \n', Rmuscle(i).nest.dof2);
    end
    if Rmuscle(i).nest.dof3 ~= []
        fprintf(fid, '                R%s \n', Rmuscle(i).nest.dof3);
    end
    if Rmuscle(i).nest.dof4 ~= []
        fprintf(fid, '                R%s \n', Rmuscle(i).nest.dof4);
    end
    fprintf(fid, '                end \n');
    fprintf(fid, '          parameters %d \n', Rmuscle(i).paranum);
    [mrow, ncolumn] = size(Rmuscle(i).para);
    dofnum = ncolumn-1;
    switch dofnum 
    case 1
        fprintf(fid, '            # i      p_i \n');
        fprintf(fid, '              %d      %12.6E \n', Rmuscle(i).para');
    case 2
        fprintf(fid, '            # i    j      p_ij \n');
        fprintf(fid, '              %d    %d       %12.6E \n', Rmuscle(i).para');
    case 3
        fprintf(fid, '            # i    j    k      p_ijk \n');
        fprintf(fid, '              %d    %d    %d      %12.6E \n', Rmuscle(i).para');
    case 4
        fprintf(fid, '            # i    j    k    l      p_ijkl \n');
        fprintf(fid, '              %d    %d    %d    %d      %12.6E \n', Rmuscle(i).para');
    end
    fprintf(fid, 'end \n\n');
end
