1 |
function [trm]=radial(ffreqs,dirs,wns,z,depth,xpos,ypos) |
---|
2 |
|
---|
3 |
%This transfer function takes the along beam radial velocities from an ADCP |
---|
4 |
%and transfers them to surface displacement for use in DIWASP. |
---|
5 |
%This was added to DIWASP on 6/22/07 |
---|
6 |
% Modified from Hoshimoto (1997) |
---|
7 |
|
---|
8 |
%Edits were made on 3/10/09 to fix a bug in the cutoff limits in the KZ |
---|
9 |
%and KZi parts of the transfer function. NOTE: there is a cutoff set of .4 |
---|
10 |
%for both the real and imaginary parts, but this ONLY applies to high |
---|
11 |
%frequencies. A cutoff is in place for the imaginary KZi at .05Hz to limit |
---|
12 |
%the effects of low frequency noise found in some data sets (at .01-.05Hz). |
---|
13 |
%This cutoff is set at a designated frequency to avoid losing variance at |
---|
14 |
%frequencies above .05Hz (20 second waves). |
---|
15 |
|
---|
16 |
|
---|
17 |
|
---|
18 |
transH=.4; % transducer face height off the bottom |
---|
19 |
dist=(xpos^2+ypos^2)^.5; %horizontal distance from origin |
---|
20 |
a=atan(dist/(z-transH));%beam angle from the vertical |
---|
21 |
|
---|
22 |
% compute the axis angle |
---|
23 |
if xpos > 0 && ypos > 0 % 1st quadrant |
---|
24 |
B= acos(xpos/dist); |
---|
25 |
elseif xpos < 0 && ypos < 0 % 3 quadrant |
---|
26 |
B= pi + acos(abs(xpos/dist)); |
---|
27 |
elseif xpos < 0 && ypos > 0 % 2nd quadrant |
---|
28 |
B= pi/2 + asin(abs(xpos/dist)); |
---|
29 |
elseif xpos > 0 && ypos < 0 % 4th quadrant |
---|
30 |
B = 3*pi/2 + asin(abs(xpos/dist)); |
---|
31 |
end |
---|
32 |
|
---|
33 |
%note: Dirs and B are both measured from the x-axis with positive angles |
---|
34 |
%being counterclockwise |
---|
35 |
|
---|
36 |
A1=(ones(1,length(dirs))*cos(a));%a vector of length dirs of the cos(a) |
---|
37 |
|
---|
38 |
|
---|
39 |
%what is set for the cutoff value to be used at high frequencies |
---|
40 |
cutoff=.4; |
---|
41 |
|
---|
42 |
%THE REAL PART |
---|
43 |
KZ=cosh(wns*z)./sinh(wns*depth); |
---|
44 |
FLKZ=ffreqs.*KZ; |
---|
45 |
|
---|
46 |
%find the cutoff and replace the values below it |
---|
47 |
cutindex=find(FLKZ < cutoff); |
---|
48 |
FLKZ(cutindex)=cutoff; |
---|
49 |
|
---|
50 |
|
---|
51 |
%THE IMAGINARY PART |
---|
52 |
KZi=sinh(wns*z)./sinh(wns*depth); |
---|
53 |
FLKZi=ffreqs.*KZi; |
---|
54 |
|
---|
55 |
%compute the max index |
---|
56 |
[maxvalue,maxindex]=max(FLKZi); |
---|
57 |
%find the cuttoff and replace the values for high frequencies |
---|
58 |
cutindex=find(FLKZi(maxindex:end) < cutoff); |
---|
59 |
cutindex=cutindex+(maxindex-1); |
---|
60 |
FLKZi(cutindex)=cutoff; |
---|
61 |
%now insert cutoff values at frequencies < .06Hz |
---|
62 |
lowindex=find(ffreqs < (.06*2*pi)); |
---|
63 |
lowvalue=FLKZi(max(lowindex)+1); |
---|
64 |
FLKZi(1:max(lowindex))=lowvalue; |
---|
65 |
|
---|
66 |
|
---|
67 |
%now put them together to create the transfer matrix |
---|
68 |
trm=(FLKZ*(sin(a)*cos(dirs-B))-(1i*FLKZi)*A1); |
---|
69 |
|
---|
70 |
|
---|
71 |
|
---|
72 |
|
---|
73 |
|
---|
74 |
|
---|
75 |
|
---|
76 |
|
---|
77 |
|
---|
78 |
|
---|