-
Notifications
You must be signed in to change notification settings - Fork 6
/
Channel_Impulse_response.m
109 lines (105 loc) · 2.31 KB
/
Channel_Impulse_response.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
tic
%% Parameters
%mu=1.325
speed=2.26e+8;
lamda=532*1e-9;
beam_width=3e-3;
source_div=20*pi/180; %radians (20o)
g=0.924; % Henyey-Greenstein constant
Roulette_threshold=1e-4;
alpha=10;
%----Harbor---------%
%{
a=0.295;
b=1.875;
c=a+b;
albedo=b/c;
%}
%-------Pure Sea----------%
%{
a=0.053;
b=0.003;
c=a+b;
albedo=b/c;
%}
%----------- Clear Ocean-----------------%
a=0.069;
b=0.08;
c=a+b;
albedo=b/c;
%}
D=20e-2; % Lens diameter
Z=20; % Distance between source and receiver(m) along Z-axis.
N=1e6;
%storeWT=zeros(N,1);
%storeAP=zeros(N,1);
cross=0;
I=0;
dist=zeros(6000,1);
imp=zeros(6000,1);
for i=1:N
d=0;
pho=Photon;
pho.wt=1;
%pho.x=0;
pho.x=beam_width*rand;
pho.z=0;
%pho.y= -beam_width/2 + ((beam_width/2)-(-beam_width/2))*rand;
pho.y=beam_width*rand;
thetaO= -source_div + (source_div-(-source_div))*rand;
phiO=2*pi*rand;
pho.mux=sin(thetaO)*cos(phiO);
pho.muy=sin(thetaO)*sin(phiO);
pho.muz=cos(thetaO);
while pho.wt~=0
s=-log(rand)/c;
pho=updateDir(pho,s);
if pho.z>=Z
x_int=pho.x+((Z-pho.z)/pho.muz)*pho.mux;
y_int=pho.y+((Z-pho.z)/pho.muz)*pho.muy;
if sqrt(x_int^2+y_int^2)<=D/2
I=I+pho.wt;
cross=cross+1;
imp(cross)=pho.wt;
%calc resi
resi=sqrt((pho.x-x_int)^2 + (pho.y-y_int)^2 + (pho.z-Z)^2);
d=d+s-resi;
dist(cross)=d;
end
break;
end
d=d+s;
%}
% find theta and phi and then change DCs
theta=Hen_Green(g);
phi=2*pi*rand;
pho=updateDC(pho,theta,phi);
% Reduce weight and check for roulette threshold
pho.wt=pho.wt*albedo;
if(pho.wt<=Roulette_threshold)
pho.wt=0;
end
end
end
I=I/N;
IR=10*log10(I);
toc
%%
time=dist/speed;
time=time(1:cross);
imp=imp(1:cross);
[time2,I]=sort(time);
imp2=imp(I);