function [out,out1,out2]=miswrite(varargin)
% MISWRITE writing of the MISS files
%
% 	   Possible commands are
%
% 	   syntax: miswrite('create [ ,.in1,.in2][ ,.chp][ ,.miss][,prop] [swap]',
%                            cas,'[dbin,sbin]',swap)
%
%
%          miswrite('WriteModel',fid,model) writes a MESH
%              see doc('model') for the model data structure format
%              The MISS groups are equal to the ProId. You can use
%              femesh/feutil to adjust those. 
%              model.Elt=feutil('selelt selface',model) can also be used
%              to select the skin of a model
%
%         miswrite('WriteChp',fid,model,def) write modes in a .CHP file
%             see doc('def') for standard modal deformation format
%			.mr(optional), .kr, .cr or .freq, .damp
%			 - if .kr is diagonal, we suppose it is freq
%			 - if .cr is diagonal or a scalar it is damp
%                        In case of diagonal form, give only .freq et .damp
%                        In the non-diagonal case, don't give .freq et .damp
%
%          The following fields in the model are used
%          model.name  is used as the root name for every MISS file
%          model.pl    material properties
%                  [MatId fe_mat('m_elastic','SI',2) Rho C] acoustic fluid 
%                  [MatId fe_mac('m_miss','SI',1) RO VP VS BETA
%          model.Stack entries
%           'info','Freq',w
%           'info','Source',data =  ...
%               struct('Node',[ns*3],'dir',[ns*3],'name','CurrentName')
%           'info','Control',data=struct('Node',[ns*3],'name','Name')
%              
% 	   Example: miswrite('create prop',cas,'dbin',1)
%
% 	   create File.IN  sample acoustic job 
% 	   		Writing of the file MISS.IN 
% 	                from the input parameters .name, (.GROUP), (.domain), 
%		        .source, .control, .w, .pl, (.post)
% 	   create prop
% 	   		Writing of BINARY files .0x.MASS , .0x.AMOR, .0x.RIGI
%			x are group numbers as specified in miss3d
%			see manual pages of miss3d for information
%
%  	   Structure of variable cas contains the fields:
%     	   	.post     : (optional) postprocessing by miss yes or no?
%		        1 = postprocessing (pression + ddl by default) , 
%		        0 = no postprocessing
%	   sbin/dbin
%			writing of binary files 'single precision' or 
%			'double precision'
%			for example calling mechanism is:
%				 miswrite('create',cas,dbin) 
%	   swap         boolean indicator if byte swapping needed
%
%
%	   See also MISREAD

% E. Balmes, JM Leclere, Auke Ditzel, R. Cotterau, G. Kergourlay


if nargin==1 & comstr(varargin{1},'cvs')
 out='$Revision: 1.19 $  $Date: 2006/01/27 14:28:50 $'; return;
end


if nargin==0 ;  help('miswrite');  return; end
[CAM,Cam]=comstr(varargin{1},1);carg=2;

%setpref('FEMLink','MissDocRoot',
%        'http://www.mssmat.ecp.fr/structures/perso/miss/notice');

% ------------------------------------------------------------------------
% Low level write commands
if comstr(Cam,'write'); [CAM,Cam]=comstr(CAM,6);

fid=varargin{carg};carg=carg+1;
if isempty(fid); fid=1;
elseif ischar(fid); 
 [wd,fname,ext]=fileparts(fid);
 fid=fopen(fid,'wb');
 if fid<0 error(sprintf('Cannot write to %s',[fname ext]))
 else; fprintf('Writing (%s) to %s\n',CAM,[fname ext])
 end
end
fname=fopen(fid);if isequal(fname,'"stdout"');fname='MISS';end
out=fid;

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
if comstr(Cam,'control'); %WriteControl % sdtweb('$Miss/node181.html')

 model=varargin{carg};carg=carg+1;
 if isstruct(model);
  data=stack_get(model,'info','Control','getdata');
 elseif isfield(model,'source');data=model.source; 
 else;data=struct('Node',[0 0 0],'dir',[1 1 1],'name','DEFAULT CONTROL'); 
 end
 if ~isempty(strfind(Cam,'skin')); % use model nodes as control points
  r1=feutil('getnode groupall',model);
  data.Node(end+[1:size(r1,1)],:)=r1;
 end

 fprintf(fid,'CONTROL %g\n',size(data.Node,1));
 if size(data.Node,2)==7
  for j1=1:size(data.Node,1); fprintf(fid,'%g %g %g\n',[data.Node(j1,5:7)]); end
 else
  for j1=1:size(data.Node,1); fprintf(fid,'%g %g %g\n',[data.Node(j1,1:3)]); end
 end
 fprintf('Wrote %i control points\n',size(data.Node,1));

 r1=stack_get(model,'info','Control','getdata');
 r1.DOF=data.Node(:,1)+.19;
 model=stack_set(model,'info','Control',r1);

 out=model;

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'exec'); %WriteExec % sdtweb('$Miss/nodexxx.html')

model=varargin{carg};carg=carg+1;

sdtw('_nb','EXEC section must be checked');
fprintf(fid,'EXEC UGTG CHAMP IMPEDANCE RFIC .6 .6\n');
fprintf(fid,'EXEC CONTROLE UDM\n');

return;
if ~isempty(stack_get(model,'info','Control'))
 fprintf(fid,'EXEC CONTROL UI\n');
else fprintf(fid,'EXEC UI\n');
end

fprintf(fid,'EXEC UGTG CHAMP IMPEDANCE UD0 FORCE\n');
fprintf(fid,'EXEC GLOBAL\nDOMAIN 1\nEXEC DIFF UTOT TTOT\n');
if  ~isempty(stack_get(model,'info','Control'))
  fprintf(fid,'EXEC CONTROL UTOT\n');
else fprintf(fid,'EXEC UTOT\n');
end


%EXEC UGTG CHAMP IMPEDANCE RFIC .3 .3 
%EXEC CONTROLE UDM 
%FIN


% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'freq'); %WriteFreq : sdtweb('$Miss/node174.html')

model=varargin{carg};carg=carg+1;
if isstruct(model);
 w=stack_get(model,'info','Freq','getdata');
elseif isfield(model,'w');w=model.w; 
else;w=1; 
end
if isempty(w); w=1;end
if length(w)<2; w(2)=w(1);end
if (std(diff(w))/mean(diff(w))>1e-4) % exact list 
 fprintf(fid,'FREQ %i\n',length(w));
 fprintf(fid,'%g ',w); fprintf(fid,'\n');
else
 fprintf(fid,'FREQ DE %g A %g PAS %g\n',w(1),w(end),diff(w(1:2)));
end
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'header'); %WriteHeader % sdtweb('$Miss/node149.html')

model=varargin{carg};carg=carg+1;
fprintf(fid,'GENER %s \nDATA\nTITRE \n%s\n',model.name,model.name);
fprintf(fid,'MAILLAGE %s.miss\nVERIF\n',model.name);fprintf(fid,'*\n');
fprintf(fid,'CHAMP\nLIRE %s.chp\nFINC\n',model.name);fprintf(fid,'*\n');
%xxx only chp at he moment

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'post'); %WritePost % sdtweb('$Miss/nodexxx.html')

%(optional) description of the post treatement commands in MISS
model=varargin{carg};carg=carg+1;
if ~isempty(stack_get(model,'info','Post'))
  fprintf(fid,['DOMAIN 0\nPOST\nFICH ddl\nMVFD LEGENDE\n' ...
               'FREQ TOUTES\nCHARGE 1\nDDL TOUS\nFINP\n' ...
               'DOMAIN 1\nPOST\nFICH pression\nSSOL\nFREQ TOUTES\n' ...
               'CHAMP 1\nDDL 1\nPOINT TOUS\nFINP\n']);
end

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'source'); %WriteSource : sdtweb('$Miss/node261.html')

model=varargin{carg};carg=carg+1;
if isstruct(model);
 data=stack_get(model,'info','Source','getdata');
elseif isfield(model,'source');data=model.source; 
else;data=struct('Node',[0 0 0],'dir',[1 1 1],'name','DEFAULT SOURCE'); 
end
if ~isfield(data,'name'); data.name='Default source name';end
fprintf(fid,'* %s\nINCI %i',data.name,size(data.Node,1));

for j1=1:size(data.Node,1);
  fprintf(fid,'SOURCE %g %g %g\n %g %g %g\n',data.Node(j1,:),data.dir(j1,:));
end
fprintf(fid,'EXEC INCIDENT\n');

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'mat'); %WriteMat : write materials

 pl=varargin{carg};carg=carg+1;
 for j1=1:size(pl,1)
  [typ,i1,i2]=fe_mat('typem',pl(j1,2));

  switch typ
  case 'm_elastic'
   switch i2
   case 1 % isotropic
     fprintf('*xxx elastic material not printed\n');
   case 2 % acoustic fluid
     fprintf(fid,'FLUI RO %g CELERITE %g\n',pl(j1,3),pl(j1,4));
   otherwise; error('Not supported m_elastic subtype');
   end  
  case 'm_miss'
   switch i2
   case 1 % RO VP VS BETA
        fprintf(fid,'MATE RO %g VP %g VS %g BETA %g\n', ...
          pl(j1,3),pl(j1,4),pl(j1,5), pl(j1,6));
     fprintf(fid,'*xxx elastic material not printed\n');
   otherwise; error('Not supported m_miss subtype');
   end  

  otherwise; error('Not supported m_* function');
  end
 end

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'model'); %WriteModel  %  sdtweb('$Miss/node178.html')

model=varargin{carg};carg=carg+1;

if ~isequal(model.Node,[1:size(model.Node,1)]);
 model.Node=feutil('getnode groupall',model);
 model=feutil('renumber',model);
 sdtw('_nb','WriteModel renumbered non-sequential nodes');
end
r1=stack_get(model,'info','OrientationNodes','getdata');
if ~isempty(r1)
 [EGroup,nGroup]=getegroup(model.Elt);
 model.Elt=feutil(sprintf('orient %s n %s %s',sprintf('%i ',1:nGroup), ...
  sprintf('%.15g ',r1.Node),r1.neg),model);
end
elt=model.Elt;node=model.Node(:,5:7);
if ~isfield(model,'name'); model.name=fname;end
if size(node,2)==7; node=node(:,5:7);end

fprintf(fid,'%s \n',model.name);
fprintf(fid,' %i %i \n',size(node,1), ...  %NumNode NumElt
              size(elt,1)-length(find(~finite(elt(:,1)))));
fprintf(fid,'(3e22.14)\n');
fprintf(fid,'%22.14e%22.14e%22.14e\n',node');
% Format conversion
  Conv={'quad4',[1 3 5 7]
        'tria3',[1 3 5]
        'mass1',1
        'hexa8',[1 3 5 7 13 15 17 19]};

[EGroup,nGroup]=getegroup(elt); i4=[];
for jGroup=1:nGroup
   [ElemF,i1,ElemP]= feutil('getelemf',elt(EGroup(jGroup),:),jGroup);
   i1=strmatch(ElemP,Conv(:,1));
   if isempty(i1); error(sprintf('%s not supported',ElemP));
   else
    cEGI = EGroup(jGroup)+1:EGroup(jGroup+1)-1;
   
    i2=zeros(length(cEGI),21);
    i2(:,Conv{i1,2})=elt(cEGI,1:length(Conv{i1,2}));
    i3=feval(ElemP,'prop');
    if i3(2); i2(:,21)=elt(cEGI,i3(2));  % Group=proid
    else;     i2(:,21)=jGroup;
    end
    i4(end+1)=i2(end);
    fprintf(fid,['%i %i %i %i %i %i %i %i %i %i %i %i %i %i %i %i %i %i %i %i' ...
      ' GR %i\n'],i2');
   end
end % jGroup

out=stack_set(model,'info','MISS_SDOM',unique(i4));
fclose(fid);

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'chp'); %WriteChp % sdtweb('$Miss/node250.html')

model=varargin{carg};carg=carg+1;
if carg<=nargin; def=varargin{carg};carg=carg+1;
else; def=model.def;
end
node=feutil('getnode groupall',model); % used nodes only
if ~isfield(model,'name'); model.name=fname;end

group=stack_get(model,'info','Group','getdata');
if isempty(group); fprintf(fid,'GROUP 1\n');
else ;             fprintf(fid,'GROUP %s\n',sprintf(' %i',group));
end 

i1=[fe_c(def.DOF,node(:,1),'ind',2); fe_c(def.DOF,[4 5 6 10:99]'/100,'ind')];
i1=unique(i1);
if ~isempty(i1); def.def(i1,:)=[];def.DOF(i1)=[];end
r1=feutil('dof2mode',def.def,def.DOF);

fprintf(fid,'MODE %i\n',size(def.def,2)); % allow for other DEFs

for j1=1:size(r1,3);
 fprintf(fid,'%i %20.12e %20.12e %20.12e\n',r1(:,:,j1)');
 fprintf(fid,'FIN\n');
end
fclose(fid);
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'imp'); %WriteImp - - - - - - - - - - - - - - - - - -

error('impedance file not implemented')

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
else; sdtw('_nb','''Write%s'' unknown',CAM);
end
% ------------------------------------------------------------------------
% Job starting commands
elseif comstr(Cam,'job'); [CAM,Cam]=comstr(CAM,4);

if comstr(Cam,'cpto');
 keyboard; % eval(sprintf('nas2up(''job%s'',varargin{carg:end})',CAM));
elseif comstr(Cam,'run');
 keyboard
end

% ------------------------------------------------------------------------
elseif comstr(Cam,'create') ; [CAM,Cam]=comstr(CAM,7);

if nargin<2 ; error('model not defined'); end
model=varargin{carg}; carg=carg+1;
if carg<=nargin; model.def=varargin{carg}; carg=carg+1; end
RunOpt=struct('OpenPerm','wb','Skin',0);

i1=strfind(Cam,'-run');if ~isempty(i1); 
 RunOpt.Run=1;CAM(i1+[0:3])='';[CAM,Cam]=comstr(CAM,1);
end
i1=strfind(Cam,'-skin');if ~isempty(i1); 
 RunOpt.Skin=1;CAM(i1+[0:4])='';[CAM,Cam]=comstr(CAM,1);
end

if isempty(CAM); CAM='MISS.IN'; end % default
[wd,fname,ext]=fileparts(CAM);if isempty(wd);wd=pwd;end
if ~isfield(model,'name'); model.name='';end
if isempty(model.name); model.name='DEFAULT';end
FileStack={fullfile(wd,[fname,ext])};
fid=fopen(FileStack{end},RunOpt.OpenPerm);

if (fid>0); fprintf('Writing %s\n',[fname ext]);
else; error('Unabled to write MISS.IN');
end

% Create FNAME.in file - - - - - - - - - - - - - - - - - - - - - - - -
if comstr(lower(ext),'.in')

   % write header, mesh  (.miss) to read and CHAMP (.chp) to read 
   miswrite('WriteHeader',fid,model); 
   miswrite('WriteFreq',fid,model);

   FileStack{end+1}=fullfile(wd,[model.name,'.miss']);
   model=miswrite('WriteModel',FileStack{end},model);
   if isfield(model,'def')&isfield(model.def,'DOF')&~isempty(model.def.DOF)
     FileStack{end+1}=fullfile(wd,[model.name,'.chp']);
     miswrite('WriteChp',FileStack{end},model); 
   else
     sdtw('_nb',sprintf('You have to create modes in %s.chp',model.name));
   end
   g=stack_get(model,'info','MISS_SDOM','getdata');
   fprintf(fid,'%s %g %s %g \n',['SDOM'],g,['GROUP'], g ); %xxx check # group

   miswrite('WriteMat',fid,model.pl);
   if RunOpt.Skin; 
     model=miswrite('writeControlSkin',fid,model);
     r1=stack_get(model,'info','Control','getdata');
     r1=feutil('addtest',r1,model);
     model=stack_set(model,'info','Control',r1);
   else; model=miswrite('writeControl',fid,model);
   end

   %miswrite('WriteSource',fid,model); % xxx check place

   %fprintf(fid,'%s\n',['KCM']); % xxx?
   fprintf(fid,'%s\n',['FINS']);

   fprintf(fid,'INTEGR TRIA 12 12 RECT 6 8\nFIND\n');
   fprintf(fid,'*\n');

   %fprintf(fid,'DOMAIN 1\nEXTE\nLIRE %s.imp\nFINE\n',model.name);
   fprintf(fid,'DOMAIN 1\n');
   miswrite('WriteExec',fid,model); %xxx
   miswrite('WritePost',fid,model); %xxx
   fprintf(fid,'FIN\n');

   out=model;FileStack{end+1}=fullfile(wd,[model.name,'_final.mat']);
   RunOpt.saved_model=FileStack{end};
   if fid~=1; fclose(fid);end % end of writing  MISS.IN

   [st,RunOpt.RelWd]=fileparts(wd);
   for j1=1:length(FileStack)
     [st,fname,ext]=fileparts(FileStack{j1});
     FileStack{j1}=[fname,ext];
   end
   RunOpt.FileStack=FileStack;
   save(RunOpt.saved_model,'model','RunOpt');
   fprintf('Final SDT/Miss model saved in : %s \n',RunOpt.saved_model);
   if isfield(RunOpt,'Run')&(RunOpt.Run)
    miss('run',wd1,FileStack);
   end
   out1=RunOpt;

% end of writing .imp and .chp
else
 error('Not a valid extension')
end

% end of writing  .miss file

  % writePROPfiles(NAME,model.struct.mass,model.struct.T, ...
  %                model.struct.freq,model.struct.damp,bintype,swap);

% ------------------------------------------------------------------------
elseif comstr(Cam,'post'); [CAM,Cam]=comstr(CAM,5);

% ------------------------------------------------------------------------
else; sdtw('_nb','''%s'' is unknown',CAM);
end % choice of RunOpt mode

% ------------------------------------------------------------------------
function writePROPfiles(BASEOFFILENAME,massmat,phi,wj,ksi,bintype,swap);
% This routine writes the .MASS .RIGI and .AMOR files needed by miss3d
% The files are binary. So far we have only treated the case of writing 
% in foreign data type. Soon, the other data types will be treated.
% Auke DITZEL, 19/06/2001 Ecole Centrale de Paris.
% size of the input files is as follows:
if (bintype=='sbin') 
  bytesize  = 8; minreclength = 8; precision = 32;
elseif (bintype=='dbin') 
  bytesize  = 8; minreclength = 16; precision = 64;
end
writeprecision  =       sprintf('int%d',precision);
sprecision      =       sprintf('float%d',precision);

if swap==1
  platform=computer;
  if comstr(platform,'LNX86')
    if (bintype=='sbin') foreign='b' ; %
    else                 foreign='s';
    end
  end
  if comstr(platform,'SGI64')|comstr(platform,'SGI')
    if (bintype=='sbin')  foreign='l';  %
    else                  foreign='a';  %
    end
  end
else  foreign = 'n' ; % thus foreign is native
end
%
mass 	= transpose(phi)*massmat*phi; mass	= diag(diag(mass));
% rigi 	= transpose(phi)*Kmat*phi; 
% or differently
% the next step gives the same result
rigi	= diag(diag(((wj*wj')*transpose(phi)*massmat*phi)));
amor	= diag(2*ksi*transpose(wj)*mass);
% write MASS file
i1	= [size(mass,1),size(mass,1),1,1,1];
zc_real = reshape(real(mass),1,i1(1)*i1(2));
zc_imag = reshape(imag(mass),1,i1(1)*i1(2));
zc	= [zc_real' zc_imag']';
zc      = reshape(zc,1,2*i1(1)*i1(2)*i1(3)*i1(4)*i1(5));

FILENAME=[BASEOFFILENAME '.MASS'];
fid = fopen(FILENAME,'w',foreign);
fseek(fid,0,-1);      
fwrite(fid,i1,writeprecision);

if (minreclength < i1(1)) 
  fwrite(fid,zeros(bytesize*(2*i1(1)-5),1),writeprecision);
  fseek(fid,i1(1)*bytesize*2,-1);
  fwrite(fid,zc,sprecision);
else
  fwrite(fid,zeros(bytesize*(2*minreclength-5),1),writeprecision);
  for ii = 1:i1(1),
    fseek(fid,ii*minreclength*precision/bytesize*2,-1);
    fwrite(fid,zc(2*i1(2)*i1(3)*i1(4)*i1(5)),sprecision);
    fwrite(fid,zeros(bytesize*2*(minreclength-i1(1))*i1(2)*i1(3)*i1(4)*i1(5),1),sprecision);
  end
end

fclose(fid);

% write RIGI file
i1	= [size(rigi,1),size(rigi,1),1,1,1];
zc_real = reshape(real(rigi),1,i1(1)*i1(2));
zc_imag = reshape(imag(rigi),1,i1(1)*i1(2));
zc	= [zc_real' zc_imag']';
zc      = reshape(zc,1,2*i1(1)*i1(2)*i1(3)*i1(4)*i1(5));
FILENAME= [BASEOFFILENAME '.RIGI'];
fid 	= fopen(FILENAME,'w',foreign);
fseek(fid,0,-1); fwrite(fid,i1,writeprecision);

if (minreclength < i1(1)) 
 fwrite(fid,zeros(bytesize*(2*i1(1)-5),1),writeprecision);
 fseek(fid,i1(1)*bytesize*2,-1); fwrite(fid,zc,sprecision);
else
 fwrite(fid,zeros(bytesize*(2*minreclength-5),1),writeprecision);
 for ii = 1:i1(1),
  fseek(fid,ii*minreclength*precision/bytesize*2,-1);
  fwrite(fid,zc(2*i1(2)*i1(3)*i1(4)*i1(5)),sprecision);
  fwrite(fid,zeros(bytesize*2*(minreclength-i1(1))*i1(2)* ...
                         i1(3)*i1(4)*i1(5),1),sprecision);
 end
end
fclose(fid);
% - - - - - - - - - - -write AMOR file
i1	= [size(amor,1),size(amor,1),1,1,1];
zc_real = reshape(real(amor),1,i1(1)*i1(2));
zc_imag = reshape(imag(amor),1,i1(1)*i1(2));
zc	= [zc_real' zc_imag']';
zc      = reshape(zc,1,2*i1(1)*i1(2)*i1(3)*i1(4)*i1(5));

FILENAME= [BASEOFFILENAME '.AMOR']; fid 	= fopen(FILENAME,'w',foreign);
fseek(fid,0,-1);      fwrite(fid,i1,writeprecision);

if (minreclength < i1(1)) 
  fwrite(fid,zeros(bytesize*(2*i1(1)-5),1),writeprecision);
  fseek(fid,i1(1)*bytesize*2,-1);
  fwrite(fid,zc,sprecision);
else
  fwrite(fid,zeros(bytesize*(2*minreclength-5),1),writeprecision);
  for ii = 1:i1(1),
    fseek(fid,ii*minreclength*precision/bytesize*2,-1);
    fwrite(fid,zc(2*i1(2)*i1(3)*i1(4)*i1(5)),sprecision);
    fwrite(fid,zeros(bytesize*2*(minreclength-i1(1))*i1(2)*i1(3)*i1(4)*i1(5),1),sprecision);
  end
end

fclose(fid);

