% function sys = fitmag(magdata,weight,heading,osysl_g)
%
%   FITMAG fits a stable, minimum phase transfer function
%   to magnitude data, MAGDATA, with a supplied frequency
%   domain weighting function, WEIGHT. Both these should
%   be VARYING matrices, with identical INDEPENDENT
%   VARIABLE values. FITMAG uses GENPHASE to generate
%   phase data, and FITSYS to do the fit.  OSYSL_G (optional) 
%   is the FRSP of the old fit.  If given, it will be
%   displayed along with the data.
%
%   See also: FITSYS, GENPHASE, INVFREQS, MAGFIT, MUSYNFIT, and MUSYNFLP.

function sys = fitmag(d,wt,heading,osysl_g)

 if nargin < 2
   disp(['usage: sys = fitmag(magdata,weight)']);
   return
 end
 if nargin == 2
   heading = 'CURVE FITTING';
 end

 [dtype,drows,dcols,dnum] = minfo(d);
 if dtype ~= 'vary'
   error(['magdata should be a VARYING matrix']);
   return
 end

 [wtype,wrows,wcols,wnum] = minfo(wt);
 if wtype ~= 'vary'
   error(['weight should be a VARYING matrix']);
   return
 end

 code = indvcmp(d,wt);
 if code ~= 1
    error('inconsistent VARYING data in magdata and weight')
    return
 end

	clg;subplot(211); 			% begin JCD CHANGE 
	 if nargin <= 3
		vplot('liv,lm',d);
 		xlabel('  data ');
	  else
		vplot('liv,lm',d,osysl_g);
 		xlabel('  data and old fit ');
	 end	
	title([heading]);
	subplot(212);
	vplot('liv,lm',wt); 
	title('wt for fit');  			%  end JCD CHANGE
 xlabel(['NOTE APPROXIMATE ORDER NECESSARY FOR FIT.....']);
 ord = input('ENTER ORDER OF CURVE FIT      ');
 while any([ isempty(ord)  (ord<0)  (floor(ord) ~= ceil(ord)) ])
     ord = input('try again - a nonnegative integer..    ');
 end

 resp = genphase(d);

 while ord >= 0
   sysh = fitsys(resp,ord,wt,1);
   sysh_g = frsp(sysh,d);
   rat_data = sbs(d,sysh_g);
   hhead = [',   W/ORDER = ' int2str(ord)];

%   vplot('liv,lm',rat_data)  %  replaced by the following code
	clg;subplot(211); 			% begin JCD CHANGE
	 if nargin <= 3
%		vplot('liv,lm',sysh_g,resp);
% 		xlabel('  1) fit    2) data ');
		vplot('liv,lm',resp,sysh_g);
 		xlabel('  1) data    2) newfit ');
	  else
%		vplot('liv,lm',sysh_g,resp,osysl_g);
% 		xlabel('  1) fit    2) data   3) oldfit ');
		vplot('liv,lm',resp,osysl_g,sysh_g);
 		xlabel('  1) data    2) oldfit   3) newfit ');
	 end
	title([heading hhead]);
	subplot(212);
	vplot('liv,lm',wt); 
	title('weight for fit');			% end JCD CHANGE

   xlabel('ENTER NEW ORDER, or NEGATIVE NUMBER TO STOP ');
   pause;
   ord = input('ENTER NEW ORDER, NEGATIVE TO STOP    ');
   while any([ isempty(ord)  (floor(ord) ~= ceil(ord)) ])
       ord = input('try again - an integer, negative to stop..    ');
   end %while any
 end %while ord
sys = sysh;
subplot(111);
%
% Copyright MUSYN INC 1991,  All Rights Reserved
