%MKCSF  Makes C code for S-function blocks using Legacy Tool
% MKCSF is a collection of calls to DO(S,CSFNAME,CSFOUT,FF,REP) with S
% smeq|smkl for MEQ|MKL S-functions, CSFNAME as '?<cfunctionname>csf',
% CSFOUT the Legacy Tool coded call to cfunctionname, FF is 's'imple or
% 'd'ouble or 'sd' replacing '?', REP as {S1 R1 S2 R2 ...} replacements in
% CSFOUT.
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
function mkcsf

 assert(~verLessThan('matlab','9.6.0'),'MEQ Simulink models now require MATLAB version>=9.6.0')
 
 smeq = legacy_code('initialize');
 smeq.IncPaths       = { pwd };
 smeq.HeaderFiles    = { 'meq.h' };
 smeq.LibPaths       = { pwd };
 smeq.HostLibFiles   = { 'meqld(''host'')' }; % will select the appropriate loader script
 smeq.TargetLibFiles = { 'meqld(''target'')' };
 
 do('init')
 fs =  's'; % build single
 fd =  'd'; % build double
 fa = 'sd'; % build all
 
 % gszr(Fx,Fb,Jy,cx,cq,cr,cs,ci,co,pp,nzx,nrx,dz)
 outfcn = {'void ?gszr('
  'FLT y1[(NZY+2)*(NRY+2)],' % Fx[(nzy+2)*(nry+2)]
  'FLT u2[2*(NRY+NZY)+4],'   % Fb[2*(nry+nzy)+4]
  'FLT u1[NZY*NRY],'         % Jy[nzy*nry]
  'FLT p1[],'                % cx[nry]
  'FLT p2[][],'              % cq[nry][nzy]
  'FLT p3[][],'              % cr[nry][nzy]
  'FLT p4[][],'              % cs[nry][nzy]
  'FLT p5,'                  % ci
  'FLT p6,'                  % co
  'FLT work1[NRY],'          % pp[nry]
  'INT (NZY+2),'             % nzx
  'INT (NRY+2),'             % nrx
  'FLT u3)'                  % dz
  };
 do(smeq,'?gszrcsf',outfcn,fs,{'NRY' 'size(p2,1)'
                               'NZY' 'size(p2,2)'});
  
 % nfdb(dnFc,Fx,nzx,nrx)
 outfcn = {'void ?nfdb('
  'FLT y1[(NZX+NRX-4)*2],' % dnFc[(nzx+nrx-4)*2]
  'FLT u1[NZX*NRX],'       % Fx[nzx*nrx]
  'INT p1,'                % nzx
  'INT p2)'                % nrx
  };
 do(smeq,'?nfdbcsf',outfcn,fs,{'NZX' 'p1'
                               'NRX' 'p2'});
  
 % asxy(zA,rA,FA,dz2FA,dr2FA,drzFA,dimA,nA,zX,rX,FX,dz2FX,dr2FX,drzFX,dimX,nX,Fx,FN,zx,rx,dc,nz,nr,dz,dr,idz,idr,Oasx,work)
 outfcn = {'INT y15 = ?asxy('
  'FLT y2[NA],'        % zA[na]
  'FLT y1[NA],'        % rA[na]
  'FLT y3[NA],'        % FA[na]
  'FLT y5[NA],'        % dz2FA[na]
  'FLT y4[NA],'        % drzFA[na]
  'FLT y6[NA],'        % dr2FA[na]
  'INT (NA),'          % na
  'INT y7[1],'         % nA
  'FLT y9[NX],'        % zX[nx]
  'FLT y8[NX],'        % rX[nx]
  'FLT y10[NX],'       % FX[nx]
  'FLT y12[NX],'       % dz2FX[nx]
  'FLT y11[NX],'       % drzFx[nx]
  'FLT y13[NX],'       % dr2FX[nx]
  'INT (NX),'          % nx
  'INT y14[1],'        % nX
  'FLT u1[NZX*NRX],'   % Fx[nzx*nrx]
  'FLT u2,'            % FN
  'FLT p2[],'          % zx[nzx]
  'FLT p1[],'          % rx[nrx]
  'INT p8,'            % dc
  'INT NZX,'           % nzx
  'INT NRX,'           % nrx
  'FLT p4,'            % dz
  'FLT p3,'            % dr
  'FLT p6,'            % idz
  'FLT p5,'            % idr
  'BOOL p7[][],'       % Oasx[nzx*nrx]
  'INT work1[3*NX])'}; % work[3*nx]
 do(smeq,'?asxycsf',outfcn,fs,{'NA'  '1'
                               'NX'  '20'
                               'NRX' 'size(p1,1)'
                               'NZX' 'size(p2,1)'});
 
 % bavx(l,a,rX,zX,vrX,vzX,nX,r,z,n)
 outfcn = {'void ?bavx('
  'BOOL Y1[N],' % l[n]
  'FLT p1,'     % a
  'FLT u1[],'   % rX 
  'FLT u2[],'   % zX
  'FLT u3[],'   % vrX 
  'FLT u4[],'   % vzX
  'INT u5,'     % nX
  'FLT u6[],'   % r[n]
  'FLT u7[],'   % z[n]
  'INT N)'};    % n
 do(smeq,'?bavxcsf',outfcn,fs,{'N' 'size(u6,1)'});
 
 % bavx2(Ol,a,rX,zX,vrX,vzX,nX,r,z,nr,nz)
 outfcn = {'void ?bavx2('
  'BOOL Y1[NZ][NR],' % Ol[nz][nr]
  'FLT p1,'          % a
  'FLT u1[],'        % rX 
  'FLT u2[],'        % zX
  'FLT u3[],'        % vrX 
  'FLT u4[],'        % vzX
  'INT u5,'          % nX
  'FLT u6[],'        % r[nr]
  'FLT u7[],'        % z[nz]
  'INT NR,'          % nr
  'INT NZ)'};        % nz
 do(smeq,'?bavx2csf',outfcn,fs,{'NR' 'size(u6,1)'
                                'NZ' 'size(u7,1)'});
 
 % fbnd(FB,rB,zB,lB,lX,Fl,rl,zl,FX,rX,zX,nX,FN,nl)
 outfcn = {'void ?fbnd('
  'FLT y1[1],'  % FB
  'FLT y2[1],'  % rB
  'FLT y3[1],'  % zB
  'BOOL y4[1],' % lB
  'BOOL y5[1],' % lX
  'FLT u1[],'   % Fl[nl]
  'FLT p1[],'   % rl[nl]
  'FLT p2[],'   % zl[nl]
  'FLT u2[],'   % FX[nx]
  'FLT u3[],'   % rX[nx]
  'FLT u4[],'   % zX[nx]
  'INT u5,'     % nX
  'FLT u6,'     % FN
  'INT NL)'};   % nl
 do(smeq,'?fbndcsf',outfcn,fs,{'NL' 'size(p1,1)'});
 
 % fl4p(Fl,Fx,kl,cl,FN,rl,zl,rA,zA,nl,nz)
 outfcn = {'void ?fl4p('
  'FLT y1[NL],' % Fl[nl]
  'FLT u1[],'   % Fx[nz*nr]
  'INT p1[],'   % kl[nl]
  'FLT p2[][],' % cl[4][nl]
  'FLT u2,'     % FN
  'FLT u3[],'   % rl[nl]
  'FLT u4[],'   % zl[nl]
  'FLT u5,'     % rA
  'FLT u6,'     % zA
  'INT NL,'     % nl
  'INT p3)'};   % nz
 do(smeq,'?fl4pcsf',outfcn,fs,{'NL' 'size(p1,1)'});
 
 % pdom(Opy,Fx,FB,Ip,OXy,nzy,nry)
 outfcn = {'void ?pdom('
  'BOOL y1[NZY][NRY],' % Opy[nzy,nry]
  'FLT u1[],'          % Fx[(nzy+1)*(nry+1)]
  'FLT u2,'            % FB
  'FLT u3,'            % Ip
  'BOOL u4[][],'       % OXy[nzy][nry]
  'INT NZY,'           % nzy
  'INT NRY)'};         % nry
 do(smeq,'?pdomcsf',outfcn,fs,{'NZY' 'size(u4,1)'
                               'NRY' 'size(u4,2)'});
 
 % bbox(x1,y1,x2,y2,lxy,x,y,nx,ny)
 outfcn = {'void ?bbox('
  'FLT y3[1],'   % x1
  'FLT y1[1],'   % y1
  'FLT y4[1],'   % x2
  'FLT y2[1],'   % y2
  'BOOL u1[][],' % lxy[nx][ny]
  'FLT p1[],'    % x[nx]
  'FLT p2[],'    % y[ny]
  'INT NX,'      % nx
  'INT NY)'};    % ny
 do(smeq,'?bboxcsf',outfcn,fs,{'NX' 'size(p1,1)'
                               'NY' 'size(p2,1)'});
 
 % bf3p1(Tyg,Tpg,ITpg,Fx,FA,FB,Opy,ry,iry,nzy,nry,circ,err)
 outfcn = {'void ?bf3p1('
 'FLT y1[NZY*NRY][3],' % Tyg[nzy*nry][3]
 'FLT y2[1][3],'       % Tpg [1][3]
 'FLT y3[1][3],'       % ITpg [1][3]
 'FLT u1[],'           % Fx[(nzy+2)*(nry+2)]
 'FLT u3,'             % FA
 'FLT u4,'             % FB
 'INT8 u2[][],'        % Opy[nzy][nry]
 'FLT p1[],'           % ry[nry]
 'FLT p2[],'           % iry[nry]
 'INT p3,'             % nzy % must be parameters for size(y1,1)
 'INT p4,'             % nry
 'BOOL u5,'            % circ
 'INT u6)'};           % err
 do(smeq,'?bf3p1csf',outfcn,fs,{'NZY' 'p3'
                                'NRY' 'p4'});
 
 % bf3p3(aPpg,aTTpg,aPg,ahqTg,ag,FA,FB,ids)
 outfcn = {'void ?bf3p3('
  'FLT y1[3],' % aPpg[3]
  'FLT y2[3],' % aTTpg[3]
  'FLT y3[3],' % aPg[3]
  'FLT y4[3],' % ahqTg[3]
  'FLT u1[3],' % ag[3]
  'FLT u2,'    % FA
  'FLT u3,'    % FB
  'FLT p1)'};  % ids
 do(smeq,'?bf3p3csf',outfcn,fs,{});

 % bf3p6(Qqg,Xq,FA,FB,rA,irA,ids)
 outfcn = {'void ?bf3p6('
  'FLT y1[1][3],' % Qqg[1,3]
  'FLT y2[1],'    % Xq[1]
  'FLT u1,'       % FA
  'FLT u2,'       % FB
  'FLT u3,'       % irA
  'FLT u4)'};     % ids
 do(smeq,'?bf3p6csf',outfcn,fs,{});
 
  % bf3p8(Bty,Fx,FA,FB,Opy,ag,rBt,err,ids,iry,nzy,nry)
 outfcn = {'void ?bf3p8('
  'FLT y1[NZY][NRY],' % Bty[nzy][nry]
  'FLT u1[],'         % Fx[(nzy+2)*(nzy+2)]
  'FLT u2,'           % FA
  'FLT u3,'           % FB
  'INT8 u4[][],'      % Opy[nzy][nry]
  'FLT u5[3],'        % ag[3]
  'FLT u6,'           % rBt
  'INT u7,'           % err
  'FLT p1,'           % ids
  'FLT p2[],'         % iry[nry]
  'INT NZY,'          % nzy
  'INT NRY)'};        % nry
 do(smeq,'?bf3p8csf',outfcn,fs,{'NZY' 'size(u4,1)'
                                'NRY' 'size(u4,2)'});
 
% bf3i1(Tyg,Tpg,ITpg,Fx,FA,FB,Opy,ry,iry,nzy,nry,gNg,IgNg,fPg,fTg,err,nN)
 outfcn = {'void ?bf3i1('
  'FLT y1[NZY*NRY][3],' % Tyg[nzy*nry][3]
  'FLT y2[1][3],'       % Tpg[3]
  'FLT y3[1][3],'       % ITpg[3]
  'FLT u1[],'           % Fx[(nzy+2)*(nry+2)]
  'FLT u3,'             % FA
  'FLT u4,'             % FB
  'INT8 u2[][],'        % Opy[nzy][nry]
  'FLT p1[],'           % ry[nry]
  'FLT p2[],'           % iry[nry]
  'INT p3,'             % nzy
  'INT p4,'             % nry
  'FLT u5[][],'         % gNg[nN][3]
  'FLT u6[][],'         % IgNg[nN][3]
  'FLT u7[],'           % fPg[1][3]
  'FLT u8[],'           % fTg[1][3]
  'INT u9,'             % err
  'INT NN)'};           % nN
 do(smeq,'?bf3i1csf',outfcn,fs,{'NZY' 'p3'
                                'NRY' 'p4'
                                'NN' 'size(u5,1)'});
 
 % bf3i2(gQg,IgQg,FNQ,nQ,gNg,IgNg,nN)
 outfcn = {'void ?bf3i2('
  'FLT y1[NQ][NG],' % gQg[nQ][3]
  'FLT y2[NQ][NG],' % IgQg[nQ][3]
  'FLT u1[],'      % FNQ[nQ]
  'INT NQ,' % nQ
  'FLT u2[][],' % gNg[nN][3]
  'FLT u3[][],' % IgNg[nN][3]
  'INT NN)'}; % nN
 do(smeq,'?bf3i2csf',outfcn,fs,{'NQ' 'size(u1,1)'
                                'NN' 'size(u2,1)'
                                'NG' 'size(u2,2)'});

 % bf3i3(aPpg,aTTpg,aPg,ahqTg,ag,FA,FB,fPg,fTg,ids)
 outfcn = {'void ?bf3i3('
  'FLT y1[3],' % aPpg[3]
  'FLT y2[3],' % aTTpg[3]
  'FLT y3[3],' % aPg[3]
  'FLT y4[3],' % ahqTg[3]
  'FLT u1[3],' % ag[3]
  'FLT u2,'    % FA
  'FLT u3,'    % FB
  'FLT u4[3],' % fPg
  'FLT u5[3],' % fTg
  'FLT p1)'};  % ids
 do(smeq,'?bf3i3csf',outfcn,fs,{});

 % bf3i4(Tdg,Tgy,Fx,FA,FB,Opy,nzy,nry,kd,fd,nd,gNg,nN,ng)
 outfcn = {'void ?bf3i4('
  'FLT y1[ND][NG],'      % Tdg[nd][ng]
  'FLT y2[NG][NZY*NRY],' % Tgy[ng][nzy*nry]
  'FLT u1[],'            % Fx[(nzy+2)*(nry+2)]
  'FLT u2,'              % FA
  'FLT u3,'              % FB
  'INT8 u4[][],'         % Opy[nzy][nry]
  'INT NZY,'             % nzy
  'INT NRY,'             % nry
  'uint32 u6[],'         % kd[nd]
  'FLT u7[][],'          % fd[2][nd]
  'INT ND,'              % nd
  'FLT u5[][],'          % gNg[nN][ng]
  'INT NN,'              % nN
  'INT NG)'};            % nG
 do(smeq,'?bf3i4csf',outfcn,fs,{'ND'  'size(u7,1)'
                                'NG'  'size(u5,2)'
                                'NN'  'size(u5,1)'
                                'NZY' 'size(u4,1)'
                                'NRY' 'size(u4,2)'});

 % bf3i6(Qqg,Xq)
 outfcn = {'void ?bf3i6('
  'FLT y1[1][3],' % Qqg[1,3]
  'FLT y2[1])'};  % Xq[1]
 do(smeq,'?bf3i6csf',outfcn,fs,{});
 
 % bfef1(Tyg,Tpg,ITpg,Fx,FA,FB,Opy,ry,iry,nzy,nry,nP,nT,err)
 outfcn = {'void ?bfef1('
 'FLT y1[NZY*NRY][NP+NT],' % Tyg[nzy*nry][nP+nT]
 'FLT y2[1][NP+NT],'       % Tpg [1][nP+nT]
 'FLT y3[1][NP+NT],'       % ITpg [1][nP+nT]
 'FLT u1[],'               % Fx[(nzy+2)*(nry+2)]
 'FLT u3,'                 % FA
 'FLT u4,'                 % FB
 'INT8 u2[][],'            % Opy[nzy][nry]
 'FLT p1[],'               % ry[nry]
 'FLT p2[],'               % iry[nry]
 'INT p3,'                 % nzy % must be parameters for size(y1,1)
 'INT p4,'                 % nry
 'INT p5,'                 % nP
 'INT p6,'                 % nT
 'INT u5)'};               % err
 do(smeq,'?bfef1csf',outfcn,fs,{'NZY' 'p3'
                                'NRY' 'p4'
                                'NP'  'p5'
                                'NT'  'p6'});
 
 % bfef3(aPpg,aTTpg,aPg,ahqTg,ag,FA,FB,ids,nP,nT)
 outfcn = {'void ?bfef3('
  'FLT y1[NP+NT],' % aPpg[nP+nT]
  'FLT y2[NP+NT],' % aTTpg[nP+nT]
  'FLT y3[NP+NT],' % aPg[nP+nT]
  'FLT y4[NP+NT],' % ahqTg[nP+nT]
  'FLT u1[NP+NT],' % ag[nP+nT]
  'FLT u2,'        % FA
  'FLT u3,'        % FB
  'FLT p1,'        % ids
  'INT p2,'        % nP
  'INT p3)'};      % nT
 do(smeq,'?bfef3csf',outfcn,fs,{'NP'  'p2'
                                'NT'  'p3'});

 % bfef6(Qqg,Xq,FA,FB,rA,irA,ids,nP,nT)
 outfcn = {'void ?bfef6('
  'FLT y1[NQ][NP+NT],' % Qqg[nq,nP+nT]
  'FLT y2[NQ],'        % Xq[nq]
  'FLT u1,'            % FA
  'FLT u2,'            % FB
  'FLT u3,'            % rA
  'FLT u4,'            % irA
  'FLT u5,'            % ids
  'INT p1,'            % nP
  'INT p2,'            % nT
  'INT p3)'};          % nq 
 do(smeq,'?bfef6csf',outfcn,fs,{'NP'  'p1'
                                'NT'  'p2'
                                'NQ'  'p3'});
 
 % bfpr(IgNg,gNg,nN)
 outfcn = {'void ?bfpr('
  'FLT y1[NN][NG],' % IgNg[nN][ng]
  'FLT u1[][],' % gNg[nN][ng]
  'INT NN)'}; % nN
 do(smeq,'?bfprcsf',outfcn,fs,{'NN' 'size(u1,1)'
                               'NG' 'size(u1,2)'});
 
 % resp(y,a,u,b,v,m,na,nu,nb,nv)
 outfcn = {'void ?resp('
  'FLT y1[M][NU+NV],' % y
  'FLT p1[][],' % a[m,na]
  'FLT u1[][],' % u[na,nu]
  'FLT p2[][],' % b[m,nb]
  'FLT u2[],'   % v[nb]
  'INT M,' % m
  'INT NA,' % na
  'INT NU,' % nu
  'INT NB,' % nb
  'INT (NV))'};       % nv=1
 do(smeq,'?respcsf',outfcn,fs,{'M'  'p3'
                               'NA' 'size(p1,2)'
                               'NU' 'p4'
                               'NB' 'size(p2,2)'
                               'NV' '1'})
 
 % bslv(ai,Je,ai1,Je1,Je0,Aie,Aei,ni,ne,niter)
 outfcn = {'void ?bslv('
  'FLT y1[NI],' % ai[ni]
  'FLT y2[NE],' % Je[ne]
  'FLT u3[],'   % ai1[ni]
  'FLT u4[],'   % Je1[ne]
  'FLT u5[],'   % Je0[ne]
  'FLT u2[][],' % Aje[ni][ne]
  'FLT u1[][],' % Aej[ne][ni]
  'INT NI,'     % ni
  'INT NE,'     % ne
  'INT p1)'};   % niter
 do(smeq,'?bslvcsf',outfcn,fa,{'NI' 'size(u3,1)'
                               'NE' 'size(u4,1)'});
 
 % iata(y,a,na,b,nb,n)
 outfcn = {'void ?iata('
  'FLT y1[NI][NI],' % y[ni][ni]
  'FLT u1[][],'     % a[na][ni]
  'INT NA,'         % na
  'FLT u2[][],'     % b[nb][ni]
  'INT NB,'         % nb
  'INT NI)'};       % ni
 do(smeq,'?iatacsf',outfcn,fa,{'NA' 'size(u1,1)'
                               'NB' 'size(u2,1)'
                               'NI' 'size(u1,2)'});

 % rtci(a,dan,ai,ero,ezo,Fx,Fq,cr,cz,Opy,Fo,Opo,dap,nr,nz,nq,no,mo)
 outfcn = {'void ?rtci('
  'FLT y1[NO][NQ],' % a[no][nq]
  'FLT y2[1],'    % dan
  'FLT u1[][],'   % ai[no][nq]
  'FLT u2[],'     % ero[1|no]
  'FLT u3[],'     % ezo[1|no]
  'FLT u4[],'     % Fx[nz*nr]    
  'FLT u5[],'     % Fq[nq]
  'FLT p1[],'     % cr[no]
  'FLT p2[],'     % cz[no]
  'INT8 u6[],'    % Opy[ny]
  'FLT u7[],'     % Fo[1|no]
  'INT8 u8,'      % Opo
  'FLT p3,'       % dap
  'INT p4,'       % nr
  'INT p5,'       % nz
  'INT NQ,'       % nq
  'INT NO,'       % no
  'INT MO)'};     % mo
 do(smeq,'?rtcicsf',outfcn,fs,{'NO' 'size(u1,1)'
                               'NQ' 'size(u1,2)'
                               'MO' 'size(u2,1)'});

 % rtci1(y,a,c,yo,nq,no,mo)
 outfcn = {'void ?rtci1('
  'FLT y1[NO][NQ],' % y[no][nq]
  'FLT u1[][],'     % a[no][nq]
  'FLT u2[],'       % c[no]
  'FLT u3[],'       % yo[mo]
  'INT NQ,'         % nq
  'INT NO,'         % no
  'INT MO)'};       % mo
 do(smeq,'?rtci1csf',outfcn,fs,{'NO' 'size(u1,1)'
                                'NQ' 'size(u1,2)'
                                'MO' 'size(u3,1)'});

 % fsgi(Q0Q,Q1Q,Q2Q,Q3Q,Q4Q,iqQ,M1q,M2q,rq,irq,rA,FA,FB,d2nFA,lX,rB,iTQ,idoq,nq,no)
 outfcn = {'void ?fsgi('
  'FLT y1[NQ1],' % Q0Q[nq+1]
  'FLT y2[NQ1],' % Q1Q[nq+1]
  'FLT y3[NQ1],' % Q2Q[nq+1]
  'FLT y4[NQ1],' % Q3Q[nq+1]
  'FLT y5[NQ1],' % Q4Q[nq+1]
  'FLT y6[NQ1],' % iqQ[nq+1]
  'FLT u1[][],'   % M1q[no][nq]
  'FLT u2[][],'   % M2q[no][nq]
  'FLT u3[][],'   % rq[no][nq]
  'FLT u4[][],'   % irq[no][nq]
  'FLT u5,'       % rA
  'FLT u6,'       % FA
  'FLT u7,'       % FB
  'FLT u8,'       % d2nFA
  'BOOL u9,'      % lX
  'FLT u10,'      % rB
  'FLT u11[],'    % iTQ[nq+1]
  'FLT p1,'       % idoq
  'INT NQ,'       % nq
  'INT NO)'};     % no
 do(smeq,'?fsgicsf',outfcn,fs,{'NQ1' 'size(u11,1)'
                               'NO'  'size(u1,1)'
                               'NQ'  'size(u1,2)'});

%shap(rgeom,zgeom,amino,epsil,kappa,delta,deltl,deltu,rrmax,zrmax,rrmin,zrmin,rzmax,zzmax,rzmin,zzmin,rq,zq,rB,zB,dr,dz,npq,noq)
 outfcn = {'void ?shap('
  'FLT  y1[NQ],'   % rgeom
  'FLT  y2[NQ],'   % zgeom
  'FLT  y3[NQ],'   % amino
  'FLT  y4[NQ],'   % epsil
  'FLT  y5[NQ],'   % kappa
  'FLT  y6[NQ],'   % delta
  'FLT  y7[NQ],'   % deltl
  'FLT  y8[NQ],'   % deltu
  'FLT  y9[NQ],'   % rrmax
  'FLT y10[NQ],'   % zrmax
  'FLT y11[NQ],'   % rrmin
  'FLT y12[NQ],'   % zrmin
  'FLT y13[NQ],'   % rzmax
  'FLT y14[NQ],'   % zzmax
  'FLT y15[NQ],'   % rzmin
  'FLT y16[NQ],'   % zzmin
  'FLT  u1[][],'   % rq
  'FLT  u2[][],'   % zq
  'FLT  u3,'       % rB
  'FLT  u4,'       % zB
  'FLT  p1,'       % dr
  'FLT  p2,'       % dz
  'INT  p3,'       % npq
  'INT  p4)'};     % noq
 do(smeq,'?shapcsf',outfcn,fs,{'NQ' 'size(u1,2)'});
               
 % vizr(Wpol,Ft0,Fx,Opy,ry,iry,rBt,dr,dz,nr,nz)
 outfcn = {'void ?vizr('
  'FLT y1[1],'   % Wpol
  'FLT y2[1],'   % Ft0
  'FLT y3[1],'   % Vp
  'FLT u1[],'    % Fx[nz*nr]
  'BOOL u2[][],' % Opy[nz-2][nr-2]
  'FLT p1[],'    % ry[nr-2]
  'FLT p2[],'    % iry[nr-2]
  'FLT u3,'      % Brt
  'FLT p3,'      % dr
  'FLT p4,'      % dz
  'INT p5,'      % nr
  'INT p6)'};    % nz
 do(smeq,'?vizrcsf',outfcn,fs,{});

 % ipmh(Ie,Jh,st,Ahd,Aed,Ahe,Aeh,Ahh,uAhh,Yd,Ie0,Jh0,sIp,nit,tol,w,nh,ne,nd)
 outfcn = {'void ?ipmh('
  'FLT y1[NE],'    % Ie[ne]
  'FLT y2[NH],'    % Jh[nh]
  'BOOL y3[1],'    % st
  'FLT p1[][],'    % Ahd[nh][nd]
  'FLT p2[][],'    % Aed[ne][nd]
  'FLT p3[][],'    % Ahe[nh][ne]
  'FLT p4[][],'    % Aeh[ne][nh]
  'FLT p5[][],'    % Ahh[nh][nh]
  'FLT p6[],'      % uAhh[nh*(nh+1)/2]
  'FLT u1[],'      % Yd[nd]
  'FLT u2[],'      % Ie0[ne]
  'FLT u3[],'      % Jh0[nh]
  'FLT u4,'        % sIp
  'INT p7,'        % nit
  'FLT u5,'        % tol
  'FLT work1[NW],' % w[nh*(nh+1)/2+8*nh+ne]
  'INT NH,'        % nh
  'INT NE,'        % ne
  'INT ND)'};      % nd
 do(smeq,'?ipmhcsf',outfcn,fa,{ 'ND' 'size(u1,1)'
                                'NE' 'size(u2,1)'
                                'NH' 'size(u3,1)'
                                'NW' 'size(p6,1)+8*size(p1,1)+size(p2,1)'});

 % ipmj(Ie,aj,st,uAjj,Aje,Aej,Ie0,aj0,Ie_,aj_,Qcj,Xc,s,z,tol,iters,w,nj,ne,nc)
 outfcn = {'void ?ipmj('
  'FLT y1[NE],'    % Ie[ne]
  'FLT y2[NJ],'    % aj[nj]
  'BOOL y3[1],'    % st
  'FLT u1[],'      % uAjj[nj*(nj+1)/2]
  'FLT u2[][],'    % Aje[nj][ne]
  'FLT u3[][],'    % Aej[ne][nj]
  'FLT u4[],'      % Ie0[ne]
  'FLT u5[],'      % aj0[nj]
  'FLT u6[],'      % Ie_[ne]
  'FLT u7[],'      % aj_[nj]
  'FLT u8[][],'    % Qcj[nc][nj]
  'FLT u9[],'      % Xc[nc]
  'FLT u10[],'     % s[nc]
  'FLT u11[],'     % z[nc]
  'FLT p1,'        % tol
  'INT p2,'        % iters
  'FLT work1[NW],' % w[nj*(nj+1)/2+4*nj+6*nc]
  'INT NJ,'        % nj
  'INT NE,'        % ne
  'INT NC)'};      % nc
 do(smeq,'?ipmjcsf',outfcn,fa,{ 'NE' 'size(u3,1)'
                                'NJ' 'size(u2,1)'
                                'NC' 'size(u8,1)'
                                'NW' 'size(u1,1)+4*size(u2,1)+6*size(u8,1)'});


 % locQ(aS,aQ,qQ,qS,NF,nQ,nO,nS,1)
 outfcn = {'void ?locQ('
  'FLT y1[NO][NS],'    % aR[nO][nS]
  'FLT u1[][],'        % aQ[nO][nQ]
  'FLT u2[],'          % raQ[nQ]
  'FLT u3[],'          % raS[nS]
  'FLT u4,'            % NF
  'INT NQ,'            % nQ
  'INT NO,'            % nO
  'INT NS,'            % nS
  'INT NAR)'};         % naR
 do(smeq,'?locScsf',outfcn,fs,{ 'NQ'  'size(u1,2)'
                                'NO'  'size(u1,1)'
                                'NS'  'size(u3,1)'
                                'NAR' '1'});

 % locQ(aR,aQ,qQ,qR,NF,nQ,1,nR,naR)
 outfcn = {'void ?locQ('
  'FLT y1[NR][NAR],'   % aR[nR][naR]
  'FLT u1[],'          % aQ[nQ]
  'FLT u2[],'          % qQ[nQ]
  'FLT u3[],'          % qR[nR]
  'FLT u4,'            % NF
  'INT NQ,'            % nQ
  'INT NO,'            % nO
  'INT NR,'            % nR
  'INT NAR)'};         % naR
 do(smeq,'?locRcsf',outfcn,fs,{ 'NQ'  'size(u1,1)'
                                'NO'  '1'
                                'NR'  'p1'
                                'NAR' 'p2'});
                              
                              
 % minQ(aqmin,qmin,aQ,qQ,sq,NF,nQ,nmin)
 outfcn = {'void ?minQ('
  'FLT y1[NMIN],'   % aqmin[nmin]
  'FLT y2[NMIN],'   % qmin[nmin]
  'FLT u1[],'       % aQ[nQ]
  'FLT u2[],'       % qQ[nQ]
  'FLT u3,'         % sq
  'FLT u4,'         % NF
  'INT NQ,'         % nQ
  'INT p1)'};       % nmin
 do(smeq,'?minQcsf',outfcn,fs,{ 'NQ'   'size(u1,1)'
                                'NMIN' 'p1'});
 
 % qint(Fi,Bri,Bzi,Brri,Brzi,Bzri,Bzzi,rx,zx,Fx,ri,zi,inM,ni,nr,nz,n)
 outfcn = {'void ?qint('
  'FLT y1[NI],' % Fi[ni]
  'FLT y2[NI],' % Bri[ni]
  'FLT y3[NI],' % Bzi[ni]
  'FLT y4[NI],' % Brri[ni]
  'FLT y5[NI],' % Brzi[ni]
  'FLT y6[NI],' % Bzri[ni]
  'FLT y7[NI],' % Bzzi[ni]
  'FLT p1[],'   % rx[nr]
  'FLT p2[],'   % zx[nz]
  'FLT u1[NZ*NR],' % Fx[nz,nr]
  'FLT u2[],'   % ri[ni]
  'FLT u3[],'   % zi[ni]
  'FLT p3[][],' % inM[n,6]
  'INT NI,'     % ni
  'INT NR,'     % nr
  'INT NZ,'     % nz
  'INT N)'};    % n
 do(smeq,'?qintcsf',outfcn,fs,{'NI' 'size(u2,1)'
                               'NR' 'size(p1,1)'
                               'NZ' 'size(p2,1)'
                               'N'  'size(p3,1)'});
 
 % bint(Fi,Fx,k,c,ni,nz)
 outfcn = {'void ?bint('
  'FLT y1[NI],' % Fi[ni]
  'FLT u1[],'   % Fx[nz,nr]
  'INT p1[],'   % k[ni]
  'FLT p2[][],' % c[4][ni]
  'INT NI,'     % ni
  'INT p3)'};   % nz
 do(smeq,'?bintcsf',outfcn,fs,{'NI' 'size(p1,1)'});
             
 % y = wclk(p)
 outfcn = {'double y1 = ?wclk(double p1, double u1)'};
 do(smeq,'?wclkcsf',outfcn,fs,{});
 
 % asgn(y,y0,u,ny,nu)
 outfcn = {'void ?asgn('
  'FLT y1[NY],' % y[ny]
  'FLT u1[],'   % y0[ny]
  'FLT u2[],'   % u[nu]
  'INT p1[],'   % iu[nu]
  'INT NY,'     % ny
  'INT NU)'};   % nu
 do(smeq,'?asgncsf',outfcn,fs,{'NY' 'size(u1,1)'
                               'NU' 'size(u2,1)'});
 
 
 % y = u.*v vdMul(n,a,b,y)
 outfcn = 'void v?Mul(INT N, FLT u1[], FLT u2[], FLT y1[N])';
 do(smeq,'?vvuvcsf',outfcn,fs,{'N' 'size(u1,1)'})
 
 % y = u./v
 % vdDiv(n,u,v,y)
 sfname = '?vdivuvcsf';
 outfcn = 'void v?Div(INT N, FLT u1[], FLT u2[], FLT y1[N])';
 do(smeq,sfname,outfcn,fs,{'N' 'size(u1,1)'})
 
 % Y = U./V vdDiv(n,U,V,Y)
 outfcn = 'void v?Div(INT M*N, FLT u1[][], FLT u2[][], FLT y1[M][N])';
 do(smeq,'?mdivuvcsf',outfcn,fs,{'M' 'size(u1,1)' 'N' 'size(u1,2)'})
 
 % y = 1./u vdInv(n,u,y)
 outfcn = 'void v?Inv(INT N, FLT u1[], FLT y1[N])';
 do(smeq,'?vinvcsf',outfcn,fs,{'N' 'size(u1,1)'})
 
 % Y = 1./U vdInv(n,U,Y)
 outfcn = 'void v?Inv(INT M*N, FLT u1[][], FLT y1[M][N])';
 do(smeq,'?minvcsf',outfcn,fs,{'M' 'size(u1,1)' 'N' 'size(u1,2)'})
 
 % y = a*U*v gemv(order=102,transA=111,m,n,a,U,ldU,v,incv,b,y,incy)
 outfcn = {'void cblas_?gemv('
  'INT (102),'
  'INT (111),'
  'INT M,'      % m
  'INT N,'      % n
  'FLT p1,'     % a
  'FLT u1[][],' % U[M][N]
  'INT M,'      % m
  'FLT u2[],'   % v[N]
  'INT (1),'    % incv=1
  'FLT (0.0),'  % b=0
  'FLT y1[M],'  % y[M]
  'INT (1))'};  % incy=1
 do(smeq,'?mvunvcsf',outfcn,fa,{'M' 'size(u1,1)' 'N','size(u1,2)'})
 
 % y = alpha*v*U gemv(order=102,transA=112,m,n,a,U,ldU,v,incv,b,y,incy)
  outfcn = {'void cblas_?gemv('
  'INT (102),'
  'INT (112),'
  'INT M,'      % m
  'INT N,'      % n
  'FLT p1,'     % a
  'FLT u1[][],' % U[M][N]
  'INT M,'      % m
  'FLT u2[],'   % v[N]
  'INT (1),'    % incv=1
  'FLT (0.0),'  % b=0
  'FLT y1[N],'  % y[N]
  'INT (1))'};  % incy=1
 do(smeq,'?mvutvcsf',outfcn,fs,{'M' 'size(u1,1)' 'N' 'size(u1,2)'})
 
 % Y = alpha*U*V gemm(order=102,transA=111,transB=111,m,n,k,a,U,m,V,k,beta,Y,m)
 outfcn = {'void cblas_?gemm('
  'INT 102,'
  'INT 111,'
  'INT 111,'
  'INT M,'        % m
  'INT N,'        % n
  'INT K,'        % k
  'FLT p1,'       % a
  'FLT u1[][],'   % U[m,k]
  'INT M,'        % m
  'FLT u2[][],'   % V[k,n]
  'INT K,'        % k
  'FLT (0.0),'    % b=0
  'FLT y1[M][N],' % Y[m][n]
  'INT M)'};      % m
 do(smeq,'?mmunvncsf',outfcn,fa,{'M' 'size(u1,1)' 'N' 'size(u2,2)' 'K' 'size(u1,2)'})
 
 % Y = alpha*U'*V gemm(order=102,transA=112,transB=111,m,n,k,a,U,ldU,V,ldV,v,Y,ldY)
 outfcn = {'void cblas_?gemm('
  'INT (102),'
  'INT (112),'
  'INT (111),'
  'INT M,'        % m
  'INT N,'        % n
  'INT K,'        % k
  'FLT p1,'       % a
  'FLT u1[][],'   % U[k][m]
  'INT K,'        % k
  'FLT u2[][],'   % V[k][n]
  'INT K,'        % k
  'FLT (0.0),'    % b=0
  'FLT y1[M][N],' % Y[m][n]
  'INT M)'};      % m
 do(smeq,'?mmutvncsf',outfcn,fs,{'M' 'size(u1,2)' 'N' 'size(u2,2)' 'K' 'size(u1,1)'})
 
 % Y = alpha*U*V'
 % ?gemm(order=102,transA=111,transB=112,m,n,k,a,U,ldU,V,ldV,b,Y,ldY)
 outfcn = {'void cblas_?gemm('
  'INT (102),'
  'INT (111),'
  'INT (112),'
  'INT M,'        % m
  'INT N,'        % n
  'INT K,'        % k
  'FLT p1,'       % a
  'FLT u1[][],'   % U[m][k]
  'INT M,'        % m
  'FLT u2[][],'   % V[n][k]
  'INT N,'        % n
  'FLT (0.0),'    % b=0
  'FLT y1[M][N],' % Y[m][n]
  'INT M)'};      % m
 do(smeq,'?mmunvtcsf',outfcn,fa,{'M' 'size(u1,1)' 'N' 'size(u2,1)' 'K' 'size(u1,2)'})
 
 do('rtwmakecfg')
end

function do(s,sfname,outfcn,flt,rep)
 persistent sall
 
 if strcmp(s,'init')
  clear sall
  load_system liulib 
  load_system mkllib 
  return
 end
 
 if strcmp(s,'rtwmakecfg')
  legacy_code('rtwmakecfg_generate',sall)
  c = fileread('rtwmakecfg.m');
  c = regexprep(c,'''(meqld\(''\w*''\))''','$1');
  % NOTE: [AM 14.10.2020] Use name with "csf" in it to allow the
  % use of a pattern rule in the Makefile and also a parallel run
  % of make. Renaming to rtwmakecfg.m will be done by the Makefile
  delete('rtwmakecfg.m');
  fid = fopen('rtwmakecfg_csf.m','w');
  fwrite(fid,c);
  fclose(fid);
  return
 end
 
 if flt
  for k = flt
   switch k
    case 's', do(s,sfname,outfcn,'',[{'?' 's' 'FLT' 'single' 'INT8' 'int8' 'INT' 'int32' 'UINT' 'uint32' 'BOOL' 'boolean'} reshape(rep',1,[])])
    case 'd', do(s,sfname,outfcn,'',[{'?' 'd' 'FLT' 'double' 'INT8' 'int8' 'INT' 'int32' 'UINT' 'uint32' 'BOOL' 'boolean'} reshape(rep',1,[])])
   end
  end
  return
 end
 
 % do it
 if iscell(outfcn), outfcn = [outfcn{:}]; end
 % user replacements
 for k = 1:2:length(rep)
  if strcmp(rep{k},'?')
   sfname = strrep(sfname,rep{k},rep{k+1});
   outfcn = strrep(outfcn,rep{k},rep{k+1});
  else
   outfcn = regexprep(outfcn,['\<' rep{k} '\>'],rep{k+1});
  end
 end
 s.SFunctionName = sfname;
 s.OutputFcnSpec = outfcn;
 fprintf('Doing %12s used in %s\n',sfname,char(find_system({'liulib' 'mkllib'},'SearchDepth',1,'BlockType','S-Function','FunctionName',sfname)))
 legacy_code('sfcn_cmex_generate',s);
 legacy_code('sfcn_tlc_generate' ,s);
 sall = [sall s];
 movefile([sfname,'.c'],fullfile('csfc',[sfname,'.c']));
  
end
