%
%   iso_apply_random.m  ver 1.2   December 20, 2011
%
function[q,fig_num]=...
    iso_apply_random(rd_power_trans_x,rd_power_trans_y,rd_power_trans_z,power_trans_x,power_trans_y,power_trans_z,fig_num)
    q=1;
    disp(' ');
    disp(' Apply steady-state random vibration base input? ');
    disp(' 1=yes  2=no ');
    irandom = input(' ');
%
    if(irandom==1)
    while(1)
        disp(' ');
        disp(' Enter input axis ');
        disp(' 1=X  2=Y  3=Z ');
        iaxis=input(' ');
%
        if(iaxis==1)
            acc_trans=power_trans_x;
            rd_trans=rd_power_trans_x;
            string1='  X-axis Input  ';
        end
%
        if(iaxis==2)
            acc_trans=power_trans_y;
            rd_trans=rd_power_trans_y;
            string1='  Y-axis Input  ';           
        end
%
        if(iaxis==3)
            acc_trans=power_trans_z;
            rd_trans=rd_power_trans_z;
            string1='  Z-axis Input  ';              
        end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
        disp(' ');
        disp(' The input PSD must have two columns:  freq(Hz) & PSD(G^2/Hz)');
        disp(' ');
        disp(' Select file input method for input PSD ');
        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);
            fid = fopen(filename,'r');
            THM = fscanf(fid,'%g %g',[2 inf]);
            THM=THM';
        else
            THM = input(' Enter the matrix name:  ');
        end
%
        PSD_in=double(THM(:,2));
        freq_in=double(THM(:,1));
        nin=length(freq_in);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
        df=0.2;
%
        clear Q1;
        clear Q2;
        Q1 = real_mult_intlog(freq_in,PSD_in,df);
        f1=Q1(:,1);
        a1=Q1(:,2);
%
        rms_base=sqrt(sum(a1)*df);
%
        clear ff;
%        
        for(kv=2:4)
            clear ab;
            Q2 = real_mult_intlog(acc_trans(:,1),acc_trans(:,kv),df);
            f2=Q2(:,1);
            a2=Q2(:,2);
%
            df2=df/2;
%
            ijk=1;
            for(i=1:length(f1))
                for(j=1:length(f2))
                    if(abs(f1(i)-f2(j))<df2)
                        ab(ijk)=a1(i)*a2(j);
                        ff(ijk)=(f1(i)+f2(j))/2.;
                        ijk=ijk+1;
                        break;
                    end
                end
            end
%
%%            size(ab)
%
            RMS=sqrt(sum(ab)*df);
%
            if(kv==2)
                psd_x=ab;
                rms_x=RMS;
            end
%
            if(kv==3)
                psd_y=ab;
                rms_y=RMS;              
            end
%
            if(kv==4)
                psd_z=ab; 
                rms_z=RMS;    
            end  
%
        end
%
        figure(fig_num);
        fig_num=fig_num+1;
        plot(freq_in,PSD_in,ff,psd_x,ff,psd_y,ff,psd_z);
        grid on;
        xlabel(' Frequency(Hz) ');
        ylabel(' Accel(G^2/Hz) ');
        set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log');
        out1=sprintf('Accel PSD %s',string1);
        title(out1);
 %
        RM=[rms_x rms_y rms_z];
 %
        if(rms_x < max(RM)*1.0e-5)
            rms_x=0;
        end
  %
        if(rms_y < max(RM)*1.0e-5)
            rms_y=0;
        end
  %
        if(rms_z < max(RM)*1.0e-5)
            rms_z=0;
        end       
 %
        outb=sprintf('Base Input %5.2g GRMS',rms_base);
        outx=sprintf('X-axis Resp %5.2g GRMS',rms_x);
        outy=sprintf('Y-axis Resp %5.2g GRMS',rms_y);
        outz=sprintf('Z-axis Resp %5.2g GRMS',rms_z);
        legend (outb,outx,outy,outz);
 %      
        fmin=min(ff);
        fmax=max(ff);
        ymax2=max(psd_x);
        ymax3=max(psd_y);
        ymax4=max(psd_z);
        ymax=ymax2;
        if(ymax<ymax3)
            ymax=ymax3;
        end
        if(ymax<ymax4)
            ymax=ymax4;
        end
 %
        for i=-10:10
            if(ymax<10^i)
                ymax=10^i;
                break;
            end
        end
%   
        ymin=ymax*0.00001;
%
        for iv=-10:10
            if fmin==(10^iv)
                fmin=(10^iv);
                break;
            end
            if fmin>(10^iv)
            else
                fmin=10^(iv-1);
                break;
            end    
        end
%
        axis([fmin,fmax,ymin,ymax]);        
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
        disp(' ');
        disp(' Acceleration Response ');
        out1=sprintf('  X-axis:  %8.3g GRMS',rms_x);
        out2=sprintf('  Y-axis:  %8.3g GRMS',rms_y);
        out3=sprintf('  Z-axis:  %8.3g GRMS',rms_z);
        disp(out1);
        disp(out2);
        disp(out3);
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
        clear ff;
%        
        for(kv=2:4)
            clear ab;
            Q2 = real_mult_intlog(rd_trans(:,1),rd_trans(:,kv),df);
            f2=Q2(:,1);
            a2=Q2(:,2);
%
            df2=df/2;
%
            ijk=1;
            for(i=1:length(f1))
                for(j=1:length(f2))
                    if(abs(f1(i)-f2(j))<df2)
                        ab(ijk)=a1(i)*a2(j);
                        ff(ijk)=(f1(i)+f2(j))/2.;
                        ijk=ijk+1;
                        break;
                    end
                end
            end
%
%%            size(ab)
%
            RMS=sqrt(sum(ab)*df);
%
            if(kv==2)
                rd_psd_x=ab;
                rms_x=RMS;
            end
%
            if(kv==3)
                rd_psd_y=ab;
                rms_y=RMS;              
            end
%
            if(kv==4)
                rd_psd_z=ab; 
                rms_z=RMS;    
            end  
%
        end
%
        figure(fig_num);
        fig_num=fig_num+1;
        plot(ff,rd_psd_x,ff,rd_psd_y,ff,rd_psd_z);
        grid on;
        xlabel(' Frequency(Hz) ');
        ylabel(' Disp(in^2/Hz) ');
        set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log');
        out1=sprintf('Relative Displacement PSD %s',string1);
        title(out1);
 %
        RM=[rms_x rms_y rms_z];
 %
        if(rms_x < max(RM)*1.0e-5)
            rms_x=0;
        end
  %
        if(rms_y < max(RM)*1.0e-5)
            rms_y=0;
        end
  %
        if(rms_z < max(RM)*1.0e-5)
            rms_z=0;
        end       
 %
        outx=sprintf('X-axis Resp %5.2g inch RMS',rms_x);
        outy=sprintf('Y-axis Resp %5.2g inch RMS',rms_y);
        outz=sprintf('Z-axis Resp %5.2g inch RMS',rms_z);
        legend (outx,outy,outz);
 %      
        fmin=min(ff);
        fmax=max(ff);
        ymax2=max(rd_psd_x);
        ymax3=max(rd_psd_y);
        ymax4=max(rd_psd_z);
        ymax=ymax2;
        if(ymax<ymax3)
            ymax=ymax3;
        end
        if(ymax<ymax4)
            ymax=ymax4;
        end
 %
        for i=-10:10
            if(ymax<10^i)
                ymax=10^i;
                break;
            end
        end
%   
        ymin=ymax*0.00001;
%
        for iv=-10:10
            if fmin==(10^iv)
                fmin=(10^iv);
                break;
            end
            if fmin>(10^iv)
            else
                fmin=10^(iv-1);
                break;
            end    
        end
%
        axis([fmin,fmax,ymin,ymax]);        
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
        disp(' ');
        disp(' Relative Displacement Response ');
        out1=sprintf('  X-axis:  %8.3g inch RMS',rms_x);
        out2=sprintf('  Y-axis:  %8.3g inch RMS',rms_y);
        out3=sprintf('  Z-axis:  %8.3g inch RMS',rms_z);
        disp(out1);
        disp(out2);
        disp(out3);
%
        disp(' ');
        disp(' Perform another random base input case? ');
        disp(' 1=yes  2=no ');
        ian=input(' ');
        if(ian==2)
           break;
        end
    end
%
    end   