CPSbopkc
Jump to navigation
Jump to search
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
