BscalcE angQ

function valaccume = bscalcE_angQ(string,depscale,depint,infile)
%Created by Erik Stassinos 3/2012
%This function uses a robust regression algorithm to calculate
%regression line, its stats, and Ed(-0)
%historically input arguments have been(-r 1ed488  # # infile)
%The numbers for the depth window were picked to give the best regression
%for light availability in Plumes and Blumes proc.  List of params is in
%/home/data65/pb/PRR/SOFT prr_params.
% -r 1ed412 1 12
% -r 1ed443 1 12
% -r 1ed490 1 12
% -r 1ed510 1 12
% -r 1ed555 1 12
% -r 1ed665 1 10
% -r 1lu412 1 12
% -r 1lu443 1 12
% -r 1lu490 1 12
% -r 1lu510 1 12
% -r 1lu555 1 12
% -r 1lu665 1 10


%==Get the cast params===
%check if up of downcast
%get stuff from castid matrix
%count number of derived params
%channels ={'1Ed412','1Ed443','1Ed490','1Ed510','1Ed555','1Ed665','1Lu412','1Lu443','1Lu490','1Lu510','1Lu5551','1Lu665'}
%scale = [1,1,1,1,1,1,1,1,1,1,1,1];
%interval = [12,12,12,12,12,10,12,12,12,12,12,10]
%s = struct('channel',channels,'scale',scale,'interval',interval)
%s(L).channel,s(L).scale,s(L).interval

infile
string
[rowindex,time_atindex,depth_atindex] = index1(infile);
infp = fopen(infile,'r');
dpcnt = 0;


while feof(infp) ~= 1
    ln = fgets(infp);
    if strncmp(ln,'<data>',6) ==1
        break
    end
    
    if strncmp(ln,'<derived_parameters>',20) ==1
        dpcnt = dpcnt+1;
    end
end
ct =0;
frewind(infp);
while feof(infp) ~=1 %used when making new list_ac9
    tu = fgets(infp);
    ct = ct+1;%counts number lines
end


%==End of get cast params==
%========================================================================
%==Check if done before==
%return boolean holder if bscal was run before
%make sure it was binned

frewind(infp);
donebf = 'n';
binned = 'n';
while ~feof(infp)
    ln = fgets(infp);
    if (strncmp(ln,'bscalc',6)) ==1
        donebf ='y';
        break
    end
    if strncmp(ln,'bin_1.0_1depth',14) ==1
        binned = 'y';
    end
end

% % % if binned == 'n'
% % %     disp('file is not binned') commented 11/15/2013
% % %     return
% % % end
%==End check if done before==
%========================================================================
%==Find index for parameter==
%if can't be found, abort
linecounter =1;
frewind(infp);
%for i = 1:ct %check for header indexes
 while ~feof(infp)   
   checkline = fgets(infp);
   
   if strncmp(checkline,'<sampled_parameters>',20) ==1

       startparam = linecounter;
   end
   if strncmp(checkline,'<derived_parameters>',20) ==1
      blankind = linecounter;
   end
   if strncmp(checkline,'<data>',6) == 1
       dataind = linecounter;
       endparam = linecounter;
   end
   if strncmp(checkline,'<filters_used>',14) ==1
       filterind = linecounter;
   end
   linecounter = linecounter +1;
   
end   

   ext = 0;
   i=1;
    %for i = 1:ct%size(varargin,2) was for pulling multiple locations
    frewind(infp);


 %make the data file
 dfile = strcat('temp',infile);
 dfd = fopen(dfile,'W+');
 for linc = 1:ct
     lllc = fgets(infp);
     if linc > dataind
         fprintf(dfd,lllc);
     end
    if linc == filterind -1
        break
    end
 end
    frewind(infp);
    fclose(dfd);
   % return
 %end make datafile   
 %==========================   
    %master loop for all channels
 %==========================   

    for count = 1:ct
        varcheck = fgets(infp);
       vari = string; %<===========
        if strncmp('0Depth',varcheck,6) ==1 %changed 12-3-12
            dep = count - startparam;
        end
        
        if strncmp('EdZPitch',varcheck,8) ==1 %changed 1-14-14
            pitch = count - startparam;
        end
        
        if strncmp('EdZRoll',varcheck,7) ==1 %changed 1-14-14
            roll = count - startparam;
        end
        
        if strncmp('LuZDepth',varcheck,8) ==1 %changed 6-26-13
            dep = count - startparam;
        end
        
        if strncmp(varcheck,'<derived_parameters>',20) ==1;
            ext =1;
        end

        if strncmp(varcheck,vari,length(vari)) ==1
            
            pull(i) = count - startparam;
            i = i+1;
        end
    end
    %end


%==End of find index==
%========================================================================
%==Read all data==
%load depth into a sepparate variable
%get size of data field

%datafilepref = strrep(infile(1),'b','binmat'); just changed
%datafile = strcat(datafilepref,infile(2:end)); just changed 5-10-2012
datafile = dfile;
m = load(datafile,'.ASC');
depth = m(:,dep);
pitch_a = m(:,pitch);
roll_a = m(:,roll);
%check if up or downcast
if mean(depth(end-5:end)) < mean(depth(1:5))
 m = flipud(m); %upcast matrix
 depth = flipud(depth);
end


%==End read all data
%========================================================================
%==Do the meat, calculate K==
%qualifiers 
%delta needs to be less than depscale
%if delta>= abs(depint) winplus =i (the deepest value for depth, input to
%function)

start = 0;
data = m(:,pull(1));
winminus =-1;
winplus =-1;
allvals = [];
valaccume = [];
%filterout =[];

for g = 1:size(m,1) % changed from length(m) 9-18-2012
    delta = m(g,dep);
    if delta < -9.0E35
        continue
    end
    if winminus == -1
    if data(g) ~= -9.9E35 && delta >= depscale
        winminus = g;
    end
    end
    
    if winplus ==-1
    if delta >= abs(depint)
        winplus = g;
    else
        continue
    end
    end
    
    %===check if window is good
    pts = 1;
    if winminus == -1 || winplus == -1 %==== added  9-18-2012 to deal with crappy casts where entire channel is flagged ES

        disp('something wrong with the data')
        data(pts) = -9.9E35;
format SHORT %changed 1-3-2013
        allvals = [allvals; 'n = ',num2str(-9.9E35),' b0 = ', num2str(-9.9E35),' b1 = ', num2str(-9.9E35),' min = ', num2str(-9.9E35),' max = ', num2str(-9.9E35),' mean = ', num2str(-9.9E35),' stdDev = ', num2str(-9.9E35),' var = ', num2str(-9.9E35),' uncertainty = ', num2str(-9.9E35),' abdev = ', num2str(-9.9E35)];

        valaccume = [valaccume;allvals]; %changed 1-3-2013
     return %continue 1-3-2013
    end


    for z = winminus:winplus

        if data(z) == -9.9E35 || depth(z) == -9.9E35 || data(z) <= 0 || pitch_a(z) >=10 || roll_a(z) >= 10
         % if pitch_a(z) >=10 || roll_a(z) >= 10
         %     pitch_a(z)
         %     roll_a(z)
         %     pause
         % end
            continue
        else
            goodpts(pts,:) = z;
        end
        pts = pts + 1;
    end
    pts = length(goodpts);

    if pts <2
        disp('error need more pts to calculate')
    end
           x = depth(goodpts);
           y = log(data(goodpts));
           %y_L = log(data(goodpts));
           n = size(x,1);
           plot((y),-x);
           xlabel('channel');
           ylabel('depth (m)');
           abdev = 0;
           a = 0;
           b = 0;
           %winplus
           %winminus
           [a,b,abdev] = medfit(x,y,pts,a,b,abdev); %do fitting of datapoints

           
          % B = robustfit(x,y)
           %fitresult = fit((x),(y),'poly1')
            intcp = exp(a); 
            slp = exp(b);

            stdval = exp(std(y));
            minval = exp(min(y));
            maxval = exp(max(y));
            meann = exp(mean(y));
            syb = 0;
            sy2 =0;
            yht = [];
            lt2 =0;
            sx2 =0;
            for lt = 1:pts%length(y)
                syb = syb+y(lt); %sum of the y
                sy2 = sy2 + y(lt)*y(lt);%sum of y squared
                sx2 = sx2 + x(lt)*x(lt);
            end

            for lt2=1:length(y)
                yht(lt2) = a + (b*x(lt2));
            end
            lt2=0;
            xb =0;
            SumS =0;    
           for lt2 =1:pts%length(y)
                SumS = SumS + (y(lt2) - yht(lt2))^2;
                xb = xb + (x(lt2) - mean(x))^2;
            end
             SumS = SumS/(pts -2);

             Sbo = sqrt(SumS*(sx2/(pts*xb)));
             varr = exp(log(stdval)^2);
             clear SumS
             tstatistic = tval((100+95)*.005,n);
             unc = exp(tstatistic*Sbo);
            allvals = [allvals; 'n = ',num2str(n),' b0 = ', num2str(intcp),' b1 = ', num2str(slp),' min = ', num2str(minval),' max = ', num2str(maxval),' mean = ', num2str(meann),' stdDev = ', num2str(stdval),' var = ', num2str(varr),' uncertainty = ', num2str(unc),' abdev = ', num2str(abdev)];
            %allvals2 = ['n = ',num2str(n),' b0 = ', num2str(slp),' b1 = ', num2str(intcp),' min = ', num2str(minval),' max = ', num2str(maxval),' mean = ', num2str(meann),' stdDev = ', num2str(stdval),' var = ', num2str(varr),' uncertainty = ', num2str(unc),' abdev = ', num2str(abdev)]
            %filterout = [filterout; string,num2str(depscale),num2str(depint)]

            valaccume = [valaccume;allvals];

if winplus && winminus ~=-1
    clear ct
    clear slope_polyfit
    clear y
    clear x
    clear intercept_polyfit
    return
end
end


delete dfile
fclose all

%==End of do the meat
%========================================================================
%==Creat the output file done by bscalcloop



function t = tval(p,df)
%adapted from Peizer & Pratt JASA, vol63, p1416
%computes t-statistic val
%there are tables to verify

 t = 0;
 clear a

 positive = p >= 0.5;
    if p>0
        p = 1.0 - p;
    else
        p = p;
    end
  if (p <= 0.0 || df <= 0)
    t = 999;
  elseif (p == 0.5)
     t = 0.0;
        elseif (df == 1)
            t = 1.0 / tan((p + p) * 1.57079633);
        elseif (df == 2)
                t = sqrt(1.0 / ((p + p) * (1.0 - p)) - 2.0);
                else	
                ddf = df;
                a = sqrt(log(1.0 / (p * p)));
                aa = a * a;
                a = a - ((2.515517 + (0.802853 * a) + (0.010328 * aa))/(1.0 + (1.432788 * a) + (0.189269 * aa) + (0.001308 * aa * a)));
                t = ddf - 0.666666667 + 1.0 / (10.0 * ddf);
                 t = sqrt(ddf * (exp(a * a * (ddf - 0.833333333) / (t * t)) - 1.0));
  end
  if t>0 
      t=t;
  else
      t=-t;
  end



  

back to C-OPS/CERBERUS Mfiles