disp(' ')
disp(' invfourier.m   ver 1.2   January 19, 2012')
disp(' by Tom Irvine   Email: tomirvine@aol.com')
disp(' ')
disp(' This program calculates the inverse of a Fourier transform')
disp(' which is pre-loaded into Matlab.  ')
disp(' ')
disp(' The Fourier transform should have less than 10000 points')
disp(' ')
disp(' It must have a constant frequency step.')
disp(' ')
disp(' It must have three columns: freq(sec)  real(units)  imag(real) ')
disp(' ')
%
clear t;
clear y;
clear z;
clear tmx;
clear tmi;
clear n;
clear dt;
clear df;
clear sr;
clear f;
clear t_real;
clear t_imag;
clear freq;
clear ff;
clear zz;
clear mu;
clear hw;
clear ms;
clear rms;
clear THM;
%
disp(' ');
disp(' Select input format: ');
disp('  1=  Two columns: freq(Hz) complex   ');
disp('  2=Three columns: freq(Hz) real imag ');
%
ifor=input(' ');
%
disp(' ')
disp(' Select file input method ');
disp('   1=external ASCII file ');
disp('   2=file preloaded into Matlab ');
file_choice = input('');
%
if(file_choice==1)
     [filename, pathname] = uigetfile('*.*');
    filename = fullfile(pathname, filename);   
%    disp(' Enter the input filename ');
%    filename = input(' ','s');
    fid = fopen(filename,'r');
%
    if(ifor==1)
        THM = fscanf(fid,'%g %g',[2 inf]);      
    else
        THM = fscanf(fid,'%g %g %g',[3 inf]);
    end    
%   
    THM=THM';
else
    THM = input(' Enter the matrix name:  ');
end
%
f=double(THM(:,1));
%
if(ifor==1)
    y_real=real(double(THM(:,2)));
    y_imag=imag(double(THM(:,2)));    
else
    y_real=double(THM(:,2));
    y_imag=double(THM(:,3));
end
%
k = length(f)
%
if(k>10000)
    out1=sprintf('\n Warning: long calculation time required for high sample number.\n');
    disp(out1)
end
%
dt=1./max(f);
%
%  Output options
%
disp(' ')
disp(' Select output option ');
choice=input(' 1=plot only   2=plot & output text file  ' );
disp(' ')
%
if choice==2
   disp(' Enter the time history output filename ');
   filename1 = input(' ','s');
   fid1 = fopen(filename1,'w');
end
%    
iflag =0;
tpi=2.*pi;
%
disp(' ');
progressbar % Create figure and set starting time
%
for  n=0:(k-1)
  progressbar(n/k) % Update figure    
%
  clear sum_c;
  clear sum_s;
  clear ccc;
  clear sss;
%
  t_real(n+1)=0.;
  t_imag(n+1)=0.;  
%  
%  arg = linspace(0,tpi*n,k);
  arg = linspace(0,tpi*n*(k-1)/k,k);   
%%
  ccc=y_real.*cos(arg)'-y_imag.*sin(arg)';
  sss=y_real.*sin(arg)'+y_imag.*cos(arg)';
%
  sum_c=sum(ccc);
  sum_s=sum(sss);
%
  t_real(n+1)=t_real(n+1)+sum_c;
  t_imag(n+1)=t_imag(n+1)+sum_s;
%
  t(n+1)=n*dt;
%
  if choice==2
        fprintf(fid1,'%11.4e %11.4e %11.4e \n',t(n+1), t_real(n+1), t_imag(n+1) );      
  end    
%
end
%
progressbar(1);
if choice==2
    fclose(fid1);     
end
%
disp(' ');
disp(' Plot Time History ');
%
    figure(1); 
    plot(t,t_real);      
    xlabel('Time (sec)');
    ylabel('Amplitude');    
    out5 = sprintf(' Inverse Fourier Transform Real');
    title(out5);  
%
    figure(2); 
    plot(t,t_imag);      
    xlabel('Time (sec)');
    ylabel('Amplitude');    
    out6 = sprintf(' Inverse Fourier Transform Imaginary');
    title(out6);  
%
    zoom on
%%%%
%
inv_real=[t',t_real'];
inv_imag=[t',t_imag'];