CPSbopkc

function inde = cpsbopkc(string,depth_width,infile)
%bopkc determines the attenuation coefficients for spectral data.  The user
%specifies the channel and depth window over which to perform the
%regression.  Originally used with binned data, it determines the k value
%for each depth z. Regression is calculated using Kd = -dlnEd(ed,lu,etc)/dz
%looks like ln(Edz(z)) = ln(Ed0-) - Kdz
%modified 1-30-2014 added pitch and roll critery for window
%modified 2-8-2014 pts to pts2 goodpts
%modified 2-8-2014 badflag line 272 and 308
%modified 2-8-2014 retain first winminus(opens depth window a couple
%modified 2-13-2014 medfit was putting NAN's for some values.  replaced NAN
%with -9.9E35
%points)
%get windows
%test for good z fit
%declare vars for regression
%find boundary conditions for Tschebyshev test
%

%[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 bopkc was run before
%make sure it was binned

frewind(infp);
donebf = 'n';
binned = 'n';
while ~feof(infp)
    ln = fgets(infp);
    if (strncmp(ln,'bopkc',6)) ==1
        donebf ='y';
        break
    end
    if strncmp(ln,'bin_1.0_1depth',14) ==1
        binned = 'y';
    end
end
% % % if binned == 'n' commented on 11/23/2013 finally getting rid of
% % % binned headers
% % %     disp('file is not binned')
% % %     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
    
    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('LuZDepth',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(varcheck,'<derived_parameters>',20) ==1;
            ext =1;
        end

        if strncmp(varcheck,vari,length(vari)) ==1 %finds the variables index
            
            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);

%==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)

%================initialize
start = 0;
data = m(:,pull(1));
winminus =1;
winplus =-1;
allvals = [];
inde = [];
nextk =0;
g =1;
am = [];
bm= [];
abdevm = [];
badline =1;
%=================

%==========================================================
%=====window loop  =====================
%===================================================
%    if strncmp(string,'1Lu665',6) ==1
%        keyboard
%    end
while g <= length(depth) %master loop counter

    delta = m(g,dep);
    if delta < -9.0E35
        badline = badline +1;
        b = -9.9E35; %%%%% 9-10-12
        inde = [inde,b];%%%%%%%%%%%%%%%%%%%%added 9-10-12
       %dbq winminus = badline;
        g = g+1; %%%%%commented 9-10-12
       continue %%%%%commented 9-10-12
    end
    
    if start == 0 && winplus == -1; %only done first through to get start location

       % winminus =g;
        if depth(winminus) < -999
            winminus = winminus +1;
        end
        if delta >= abs(depth(winminus) + depth_width/2) %for starting out.
            nextk = g;
        end

        start = nextk; %+ depth_width/2

          if delta >= abs(depth(winminus) + depth_width) %if passing half window size.
               winplus = g;
              % disp('breaking')
               pause
               g = g+1;
               
               break
           else
              % disp('first incr')

 %=g; %playing with this 2-15-2013 wat just start:
               g =g+1;
               b = -9.9E35;
               inde = [inde,b];
               size(inde);
               continue
         end
    end

    if start ~= 0 %if not first time through the loop

        delta = depth_width/2;
        condition = depth(start) - delta; %checking for this "window min"
        ptr =start;
        
        while ptr ~= 0 %starts in center of bin and works down
          if winplus ~= -1 %retians first winminus from first go-through
            if depth(ptr) <= condition
            winminus = ptr;
            break
            end
          end
            ptr = ptr -1;
         end
        ptr =start; %resets ptr from center and works up.
        condition = depth(ptr) + delta;
        while ptr <= length(depth) %+1
            if depth(ptr) >= condition
                winplus = ptr;
                break
            end
            ptr = ptr +1;
        end
        start = start +1;
    end

%=================================================    
%=================end window loop===========
%=================================================
    %===check if window is good ========

    pts = 1;
    badflag =0;
    for z = winminus:winplus
        if data(z) == -9.9E35 || depth(z) == -9.9E35 || data(z) <= 0 || pitch_a(z) >=10 || roll_a(z) >= 10

  %          badflag =1;
  %          b = -9.9E35; %uncommented 1-15-13
            %inde = [inde,b]; %uncommented 1-15-13
           % continue
%            pitch_a(z)
%            roll_a(z)
%            depth(z)
%            data(z)
%            keyboard
        else
            goodpts(pts,:) = z;
        end
        pts = pts + 1;
    end
    pts = length(goodpts);
    if pts <2 %error if not enough points
        keyboard
        if badflag ==1
            %continue
        else
        disp('error need more pts to calculate')
        keyboard
        end
    end
    goodpts(find(goodpts ==0)) = [];
    pts2 = length(goodpts); %actual number of good points after cleanup up from line above
    x = depth(goodpts); % declare x and y for gregression
    y = log(data(goodpts));
    %y_L = log(data(goodpts));
    n = size(x,1);
   % plot((y),-x);  Hanging on this ?????????????????? 10/28/2013
   % xlabel('channel');
   % ylabel('depth (m)');
    [abdev,a,b,siga,sigb] = deal(0);

%==============run the medfit script ======================
%==========================================================
% % %if badflag ==0

    [a,b,abdev] = medfit(x,y,pts2,a,b,abdev);
    
    
    if winplus - winminus < 10%strncmp(string,'1Lu665',6) ==1 so computation window does not exceed winplus.
        %     winminus
        %     winplus
        % keyboard
        b = -9.9E35;
    end
% % %else
% % %  b = -9.9E35;
% % %end
goodpts = []; %clear this for next iteration of loop
if delta < -9.0E35
    b = -9.9E35;
    pause
end

%=====assign values gotten from medfit============
%======================
    if isnan(b) %added 2-12-2014
        b = -9.9E35;
    end
    inde = [inde,-b];

%======================

   % disp('second incrim')
    g = g+1;

end   


%pause(5)
 

back to C-OPS/CERBERUS Mfiles