Eriks at 20:10, 13 February 2014

2014-02-13T20:10:09Z

Eriks at 20:09, 13 February 2014

2014-02-13T20:09:36Z

Eriks at 20:07, 13 February 2014

2014-02-13T20:07:57Z

Eriks: Created page with "

function valaccume = bscalcE_angQ(string,depscale,depint,infile) %Created by Erik Stassinos 3/2012 %This function uses a robust regression algorithm to calculate %regress..."

2014-02-13T20:05:42Z

Created page with "<pre>function valaccume = bscalcE_angQ(string,depscale,depint,infile) %Created by Erik Stassinos 3/2012 %This function uses a robust regression algorithm to calculate %regress..."

New page

<pre>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

<\pre>

← Older revision		Revision as of 20:10, 13 February 2014
Line 371:		Line 371:


−	<\pre>	+	</pre>

	back to [[C-OPS/CERBERUS Mfiles]]		back to [[C-OPS/CERBERUS Mfiles]]

← Older revision		Revision as of 20:07, 13 February 2014
Line 371:		Line 371:

	<\pre>		<\pre>
		+
		+	back to [[C-OPS/CERBERUS Mfiles]]

BscalcE angQ - Revision history

Eriks at 20:10, 13 February 2014

Eriks at 20:09, 13 February 2014

Eriks at 20:07, 13 February 2014

Eriks: Created page with "
function valaccume = bscalcE_angQ(string,depscale,depint,infile) %Created by Erik Stassinos 3/2012 %This function uses a robust regression algorithm to calculate %regress..."

@@ Line 1: / Line 1: @@
-<pre>function valaccume = bscalcE_angQ(string,depscale,depint,infile)
+<pre>
 %Created by Erik Stassinos 3/2012
 %This function uses a robust regression algorithm to calculate
@@ Line 370: / Line 371: @@
-<\pre>
+  <\pre>
 back to [[C-OPS/CERBERUS Mfiles]]

BscalcE angQ - Revision history

Eriks at 20:10, 13 February 2014

Eriks at 20:09, 13 February 2014

Eriks at 20:07, 13 February 2014

Eriks: Created page with "function valaccume = bscalcE_angQ(string,depscale,depint,infile) %Created by Erik Stassinos 3/2012 %This function uses a robust regression algorithm to calculate %regress..."

Eriks: Created page with "
function valaccume = bscalcE_angQ(string,depscale,depint,infile) %Created by Erik Stassinos 3/2012 %This function uses a robust regression algorithm to calculate %regress..."