CPSdoKQ

function CPSdoKQ(infile,type,channel,depth_int,trend,variation)
%(infile,type,channel,depth_int,trend,variation)
%This function qualifies the incident illumination of a spectral channel
%during sampling over a given depth interval and writes a flag on a new
%data line
%The user chooses the channel and window width(m).  For each line the
%program calculates teh standard deviation and mean of the first difference
%(M1stdiff) for teh group of points in the interval z+-
%depth_window_width/2
%A new line will be placed in <derived_parameters> KQ-
%cmd example -m 2ed488 10 0.05 0.02 infile outfile
%compare the M1stdiff/mean and sigm/mean of data in window +- 5m to 0.02
%and 0.05 respectively for channel 2ed488
%1-15-2014 needed to change sigma to sigm-matlab function conflict


%===SET PARAMETERS===

infid = fopen(infile,'r');
dataname = strcat('KQCmat',infile);
dataname = strrep(dataname,'.lcd','');

%===MAKE DATA FILE FOR INPUT FILE===


frewind(infid);

datafid = fopen(dataname,'w');
while ~feof(infid)
    testline = fgets(infid);
    if strncmp(testline,'<data>',6) ==1 %start of data section
        while strncmp(testline,'<filters_used>',14) ~=1
            testline = fgets(infid);
              if strncmp(testline,'<filters_used>',14) ==1
                  break
              end
            fprintf(datafid,testline);
        end
    end
end
fclose(datafid);
frewind(infid);


%===GET CASTID MATRIX FROM FILE AND SPLIT DATA===
%fid = fopen(infile);

[in,mer,dep] =  index1(infile);

mat = load(dataname);

     clear ij
     clear downcastmatrix
     clear upcastmatrix
     clear diffdn
     clear diffup
             width = length(mat(1,:));
             %for ij = in{1}:1:in{2};
             for ij = 1:1:in{2}
             %downcastmatrix(ij-(in{1}-1),:) = mat(ij,1:width); %starts from 1 and goes to in2
              downcastmatrix(ij,:) = mat(ij,1:width);
             end
            % for r = in{3}:in{4};
            % upcastmatrix(r-(in{3}-1),:) = mat(r,1:width);
            % end
%mat = mat(in{1}:end,:); %index mat appropriatedly to match downcastmatrix
%===FIND CHANNEL INDEX===

lineind =0;
while ~feof(infid)
     line = fgets(infid);
                 if strncmp(line,'<sampled_parameters>',20) ==1
                     while strncmp(line,'<data>',6) ~= 1
                         lineind = lineind +1;
                         line = fgets(infid);
                         if strncmp(line,'<derived_parameters>',20) == 1
                              lineind = lineind -1;
                         end
                         if strncmp(channel,line,5) ==1
                         chanind = lineind;
                         end
                     end
                 end
end


%===INDEX OF DEPTH ONLY===

lineind =0;
frewind(infid)
while ~feof(infid)
     line = fgets(infid);
                 if strncmp(line,'<sampled_parameters>',20) ==1
                     while strncmp(line,'<data>',6) ~= 1
                         lineind = lineind +1;
                         line = fgets(infid);
                         if strncmp(line,'<derived_parameters>',20) == 1
                              lineind = lineind -1;
                         end
                         if strncmp('LuZDepth',line,8) ==1 %changed 12-3-12
                         depthind = lineind;
                         end
                     end
                 end
end
depthind;
%===INDEX OF TIME ONLY===

lineind =0;
frewind(infid);
while ~feof(infid)
     line = fgets(infid);
                 if strncmp(line,'<sampled_parameters>',20) ==1
                     while strncmp(line,'<data>',6) ~= 1
                         lineind = lineind +1;
                         line = fgets(infid);
                         if strncmp(line,'<derived_parameters>',20) == 1
                              lineind = lineind -1;
                         end
                         if strncmp('FrameCount',line,10) ==1 %changed 12-3-12 0prr_record
                         %if strncmp('1mer_time',line,9) ==1 %changed 12-3-12 0prr_record
                         timeind = lineind;
                         end
                     end
                 end
end
timeind;
dbstop if error


data = downcastmatrix(:,chanind); %specific channel data
time = downcastmatrix(:,timeind);
depth = downcastmatrix(:,depthind);

%===NEED TO SORT DATA TO BE MONOTONIC AND THEN RETURN WHEN OPERATION DONE==

% unsorted_ind = 1:length(downcastmatrix);%index to sort back to
% [depth_sorted,ds_ind] = sort(depth);%sorting just the depth and getting index
% sorted_DCmatrix = downcastmatrix(ds_ind,:);%sorting all according index of sorted depth
% new_unsort_ind(ds_ind) = unsorted_ind; %values at sorted index = unsorted index
% %new_unsorted = sorted_DCmatrix(new_unsort_ind,:);%should be same as original downcast matrix
% 
% %===Stuff is now sorted=== but currently not used
% 
% timesort = time(ds_ind,:);
% datasort = data(ds_ind,:);
% depthsort = depth_sorted;
% downcastmatrixsort = sorted_DCmatrix;

%===SORT WHOLE MATRIX JUST IN CASE=== not going to work with both up and
%downcast

%wholeunsort_ind = 1:length(mat);
%[depth_whole_sort,dsw_ind] = sort(mat(:,depthind));
%sorted_wholemat = mat(dsw_ind,:);
%new_unsortind(dsw_ind) = wholeunsort_ind;



%=====LOAD DATA AND DO=====
%================================================================
%DoKQ
%Condition                                                  Flag
%contains trapflag                                          14
%outside of cast index                                      13
%calculated M1stdiff > max M1stdiff and sig > sigmx        12
%calculated M1stdiff > max M1stdiff                         11
%calculated sig > sigmx                                    10
%point is less than winmin                                  1
%good data                                                  0
%================================================================
window = depth_int/2;
%bflag = {'notfound'};
%aflag = {'notfound'};
bflag =5; %==notfound
aflag =5; %==notfound
%for row =1:in{2} +1 %initialize result and sigm
%    sigm(row,:) =5;
%    result(row,:) =-5;
%end
r = 1;%in{1}; %originally r = {in} because it was working with entire depth matrix not pre-timmed to castID such as here
%index =1;
while r <= in{2}

   if data(r)== -9.9e35 || depth(r) == -9.9e35
       result(r) = -5;
       sigm(r) = -9.9E35;
       r =r+1;
       continue
   end
   %lower window(shallow)

   bcondition = depth(r) - window;

     for pt =r:-1:1%count down from top index to first point of data  %(length(downcastmatrixsort(:,depthind)-1))
       delta = depth(pt);%downcastmatrix(pt,depthind);

        if delta <= bcondition
            disp('===========================')
            disp('found a winmin')
           winmin = pt;
           bflag =1; %found
           delta_b = delta;
           bcondition;

           break
        else

            %disp('continueing')
            continue
        end
     end
        if bflag ==5;
            %result(r)=-9.9E35;
            result(r) =5;%NOTFOUND
            r = r+1;

            continue
        end

     %upper window(deeper)
     acondition = depth(r) + window;
     for pt = r:in{2}%length(downcastmatrixsort) %goes down in depth to find other window point 

         delta = mat(pt,depthind); %beacuse downcastmatrix is offset by index
         %if acondition > mat(in{2}-in{1},depthind)
          if acondition > mat(in{2}-1,depthind)
             r = r+1;
             result(r) = 5;
             continue
             disp('overbounds')
         end
         
         if delta >= acondition
             disp('found a winmax')
             winplus = pt
             aflag = 1;
             acondition;
             delta_a = delta;
             disp('==========================')
             break
         else
             continue
         end
     end
         if aflag==5
             %result(r)=-9.9E35;
             disp('aflag 5')
             result(r)=5; %NOTFOUND
             r = r+1;
             continue
         end
 
 
 %================got to here there was an "end" for the while < loop  
format long
%===DO M1SDIFF AND SIGM===
nint = winplus - winmin +1;
fsum =0;
tsum =0;
qout =0;
S =0;
    z = winmin;
    while z<winplus
        fsum = fsum + data(z+1) - data(z);
        tsum = tsum + time(z+1) - time(z);
        qout = qout + abs(fsum/tsum);
        z = z +1;
    end

for z = winmin:winplus
    ex = data(z);
    S = [S,ex];
    if z == winmin
        S(1) = ex;
    end
end

sigm(r) = std(S)./mean(S);
result(r) = ((1/nint)*qout)/mean(S);

r = r+1;
end


%sigm_unsort = sigm(new_unsort_ind);
%result_unsort = result(new_unsort_ind);
%===fill in column with appropriate flag===

%appline = zeros(size(result));
result_prev = result;
for i = 1:in{2}%1:length(result_unsort)

    if result(i) > trend && sigm(i) > variation && result(i) ~= 5.00
        appline(i) = 12;
    elseif result(i) > trend && result(i) ~= 5.00
            appline(i) = 11;
        elseif sigm(i) > variation %&& result(i) ~= 5.00
%                 if i == 56
%                     sigm(i)
%                     variation
%                     keyboard
%                 end
                appline(i) =10;
            elseif result(i) < -9.0e35%or outside of downcast
                    %keyboard
                    appline(i) =13;
                elseif result(i) == 5 %&& result(i) ~= 5.00
                        appline(i) =1;
                    elseif result(i) == -5 %&& result(i) ~= 1.00
                            appline(i) =14;
                        else
%                             if result(i) == 1.00
%                                 appline(i) = result(i);
%                             else
                                appline(i) =0;
%                            end
     end
  
end


%outside of downcast bounds.. just flag 
 for it = in{2}+1:length(mat)
     appline(it) = 13;
 end

%===MAKE DATA FILE===

QCdata = load(dataname);
y = size(QCdata,1);
x = size(QCdata,2);
QCdata_new = zeros(y,x+1);
QCdata_new(1:y,1:x) = QCdata;
QCdata_new(:,end) = appline';

%===INSERT NEW DATA INTO NEW LCD===

frewind(infid);
qclcd = fopen(strcat('K',infile),'w+');
while ~feof(infid)
       pos = fgets(infid);
       if strncmp(pos,'<data>',6) ==1 %start of data section
        break
       end
       fprintf(qclcd,pos);       
end

%===INSERT NEW DERIVED PARAMS AND DATA===
newparam = strcat('kq-',channel);
fprintf(qclcd,'%s\n',newparam);
fprintf(qclcd,'%s\n','<data>');
%===overwrite old matrix
dlmwrite(dataname,QCdata_new,'delimiter',' ','precision',6); %overwrite old matrix file

%orcmat = load('akrmP101216C_nhd.lcd');  %%%=============was testing here
%nwcmat = load('KQCmatrqcmmP101216C');
%minus = orcmat(:,end-2) - nwcmat(:,end);
%nzro = find(minus ~= 0)
%keyboard

datafile = fopen(dataname);
 while ~feof(datafile)
   %for loop = 1:length(QCdata_new)
        datln = fgets(datafile);
        fprintf(qclcd,datln);
 end


%===ADD TO FILTERS USED ON NEW LCD===

while ~feof(infid)
       lcdline = fgets(infid);
       if strncmp(lcdline,'<filters_used>',14) ==1
           fprintf(qclcd,'%s\n','<filters_used>');
           while ~feof(infid)
               filter = fgets(infid);
               fprintf(qclcd,'%s',filter);
           end
       end
 end
       %domathstring =strcat(['domath',type,channel1,channel2,num2str(number),outfieldname]);
       %fprintf(newlcdfile,'%s\n', domathstring);
       fprintf(qclcd,'%s','')
       fprintf(qclcd,'%s','doKQ ');
       fprintf(qclcd,'%s',' ');
       fprintf(qclcd,'%s',type);
       fprintf(qclcd,'%s',' ');
       fprintf(qclcd,'%s',channel);
       fprintf(qclcd,'%s',' ');
       fprintf(qclcd,'%s',num2str(depth_int));
       fprintf(qclcd,'%s',' ');
       fprintf(qclcd,'%s',num2str(trend));
       fprintf(qclcd,'%s',' ');
       fprintf(qclcd,'%s',num2str(variation));

fclose('all');

back to C-OPS/CERBERUS Mfiles