%% read in all original model equations
[n11,n12] = size(['Welf1';'Welf2';'Welf ']);
[n21,n22] = size(exovarmat);
[n31,n32] = size(endovarmat); 
nshk = size(shklist,2); 

maxindex  = find([n12,n22,n32] == max([n12,n22,n32]));
if size(maxindex,2)>1 
    maxindex = maxindex(1);
end

if maxindex == 3
    varmat = endovarmat;
    varmat(n31+1,1:n12) = 'Welf1';     
    varmat(n31+2,1:n12) = 'Welf2';    
    varmat(n31+3,1:n12) = 'Welf ';    
    varmat(n31+4:n31+3+n21,1:n22) = exovarmat(1:n21,1:n22);        
elseif maxindex == 2
    varmat = exovarmat;
    varmat(n21+1,1:n12) = 'Welf1';     
    varmat(n21+2,1:n12) = 'Welf2';    
    varmat(n21+3,1:n12) = 'Welf ';    
    varmat(n21+4:n21+3+n31,1:n32) = endovarmat(1:n31,1:n32);        
end    

equations = rawfile(modline+1:rulestartline(1));
equations = [equations; rawfile(startline:endmodline-1)];
equationsnew = equations;
counter = 0;
for iline = 1:size(equations,1)
    if isempty(strfind(equations{iline},'//'))
        counter = counter+1;
        equationsnew{counter} = equations{iline};
    end    
end
Nlines = counter;
equations = cell(Nlines,1);
for iline = 1:Nlines
    equations{iline} = equationsnew{iline};
end

%% for each model equation replace variables by steady state _ss notation

for iline = 1:size(equations,1)
    for i = 1:size(varmat,1)
        equations{iline} = strrep(equations{iline},[deblank(varmat(i,:)),'(-1)'],deblank(varmat(i,:)));
        equations{iline} = strrep(equations{iline},[deblank(varmat(i,:)),'(+1)'],deblank(varmat(i,:)));
        equations{iline} = strrep(equations{iline},[deblank(varmat(i,:)),'(1)'],deblank(varmat(i,:)));
    end;
end

counter = 0;
for iline = 1:size(equations,1)

    if ~isempty(strfind(equations{iline},'='))
        counter = counter+1;
        equations{iline} = strrep(equations{iline},'=','- (');
        equations{iline} = ['fval(',num2str(counter),') = ',equations{iline}];
        if isempty(strfind(equations{iline},';'))
            linebreak_indx = 1;
            linebreak_iline = iline;        
            while linebreak_indx == 1
                equations{linebreak_iline} = [equations{linebreak_iline},' ...'];        

                if ~isempty(strfind(equations{linebreak_iline+1},';'))
                    linebreak_indx = 0;
                else
                    linebreak_iline = linebreak_iline+1;
                end
            end    
        end
            
    end
    equations{iline} = strrep(equations{iline},';',');');

    for index = 1:nshk
        equations{iline} = strrep(equations{iline},char(shklist(index)),'0');
    end
end


%%

f2in = fopen('template_for_ramsey_endoss.m','r');
f2out = fopen([outfilename,'_endoss.m'],'w');

while ~feof(f2in)
    s = fgets(f2in);
    s = strrep(s,'RAMSEY_FILE',outfilename);
    fprintf(f2out,'%s',s);
end;

fclose(f2in);

for iline = 1:Nlines
  fprintf(f2out,'%s \n',char(equations(iline)));
end;

fclose(f2out);
