%
function [f]=m16_par(x,itr,ip)
%
%setting pathname to current directory
pathname=[pwd  '/'];
path1=pathname;
%
%%% Input parameters
%
p=[2000000.;1.573e11;9.67e10;5680.;3780.;3780.].*x;
    
namefiletowrite = ['parm.dat'];
fidw = fopen([pathname,namefiletowrite],'w');

fprintf(fidw,'PELAS*   8               %2.4e\n',p(1));
fprintf(fidw,'MAT1*    2               %2.4e                      0.25\n',p(2));
fprintf(fidw,'+       %5.1f\n',p(4));
fprintf(fidw,'MAT1*    3               %2.4e                      0.25\n',p(3));
fprintf(fidw,'+       %5.1f\n',p(5));
fprintf(fidw,'MAT1*    4               9.5600e+10                      0.25\n');
fprintf(fidw,'+       %5.1f\n',p(6));

fclose (fidw);

%nastran command
unix(['msc2006 nastran m16_pp_mode.bdf scr=yes batch=no memory=128mb']);

% ======================
% extracting frequencies
% ======================
namefiletoopen =['m16_pp_mode.f06'];

y = [pathname,namefiletoopen];
fid=fopen(y);
str='R E A L   E I G E N V A L U E S';
len_str=length(str);
found = 0;
while found==0
    tline = fgetl(fid);
    for(t=1:length(tline)-length(str)+1)
        countmatch = 0;
        for(count=1:length(str))
            if(str(count)==tline(t+count-1))
                countmatch=countmatch+1;
            end
        end
        if countmatch==len_str
            found =1;
            break,end
    end
end
 
tline = fgetl(fid);
tline = fgetl(fid);
freq = fscanf(fid, '%g',[7 16]);

%backup parameter file
unix(['cp parm.dat parm_' num2str(itr) '.dat']);
unix(['rm *.f04* *.xdb *.log *.pch parm.dat']);
 
% ================
% extracting modes 
% ================
namefiletoopen =['phi.op4'];
y = [pathname,namefiletoopen];
fid = fopen(y);
status=fseek(fid,0,'eof');EOF=ftell(fid);currentFPI=fseek(fid,0,'bof');

% reading MX size in first line
line=fgetl(fid);currentFPI=ftell(fid);
matrixdata=str2num(line(1:8*2));
PHI=zeros(fliplr(matrixdata));
while currentFPI<EOF
    % reading data position
    line=fgetl(fid);currentFPI=ftell(fid);
    mxPosition=str2num(line);
    if mxPosition(2)<=size(PHI,1) & mxPosition(1)<=size(PHI,2) % it always gives extra useless number in num_col+1
        counter=1;clear mxData
        while counter<=mxPosition(3)
            line=fgetl(fid);currentFPI=ftell(fid);
            newline=double(line);
            cc=0;
            for ijk=16:16:length(newline)-1;
                if newline(ijk+1)~=double(' ')
                    newline=[newline(1:ijk+cc) double(' ') newline(ijk+1+cc:end)];                
                    cc=cc+1;
                end
            end
            line=char(newline);
            data=str2num(line);
            if length(data)==5
                mxData(counter:counter+4)=data;
                counter=counter+5;
            else
                mxData(counter:counter+length(data)-1)=data;
                counter=counter+length(data);
            end
        end
        PHI(mxPosition(2):mxPosition(2)+mxPosition(3)-1,mxPosition(1))=mxData';
    else
        break    
    end
end
fclose (fid);
%
% selecting natural frequencies
f=[freq(5,7:16)];
%
end
