clear
clc
format long e
syms x0 x r0 r xx0 yy0  xx yy
global han muu dert jie fa a s I  bianchang_1 bianchang_2 h jie_x jie_y
aa=load('J1_47.dat','r');
han=aa;
ttt=0;
bianchang_1=1000;         % A邊 8.862279616441811e+001  長(zhǎng)邊
bianchang_2=1000;         % A邊 8.862279616441811e+001  短邊
I=-200;

muu=4*pi*1e-7;
fre=0;
%dert=input('電導(dǎo)率:');
dert=[1/30 0.1 1/500];
h=[30 10];
jie_z=input('接收距地面的距離:');
jie_x=input('接收x坐標(biāo):');
jie_y=input('接收y坐標(biāo):');
fa=input('發(fā)射距地面的距離:');

for ii=1:3
    b(ii)=sqrt(i*(2*pi*fre)*muu*dert(ii)); %第ii層波數(shù)
    d=sqrt(x0.^2+b(ii).^2); %Uj
    u(ii)=d;
end
%h=input('厚度:'); %厚度

mmm0=coth(u(1)*h(1)+acoth(u(1)/u(2)*coth(u(2)*h(2)+acoth(u(2)/u(3)*1))));

%f0=x0.*(1+exp(-(jie+fa).*x0).*(x0.*mmm0-u(1))./(x0.*mmm0+u(1)));     %現(xiàn)在線框距地面0.05m  發(fā)射和接收在同一平面

f0=-i*(2*pi*fre)*x0.*exp(x0.*(-jie_z)).*(1+(x0.*mmm0-u(1))./(x0.*mmm0+u(1)));%電性源半航空
a=-3.05078187595e+00; %常數(shù)
s=1.10599010095e-01;  %常數(shù)
sum0=0;
for n=1:47
   lannn0=a+(n-1).*s;
   lannni0=(1./r0).*power(10,lannn0); %被積變量
   ff0=subs(f0,x0,lannni0);
   z0=ff0.*han(n); %被積函數(shù)與漢克爾函數(shù)系數(shù)相乘
   sum0=sum0+z0;
end

g0=sum0./r0.^2;


r_1=sqrt((jie_x-xx0)^2+(jie_y).^2); % A邊 

g_1=subs(g0,r0,r_1); % A邊的積分函數(shù)



%AA的simpsion積分
a_0=-bianchang_1;%input('a_1='); %邊長(zhǎng)
b_0=bianchang_1;%input('b_1='); %邊長(zhǎng)
snn2_1=IntGauss('g_1',a_0,b_0,7)
AA0=(jie_y).*snn2_1 


%I=input('I=');

Hz0=I.*(AA0)./(4*pi)
Hz_0=abs(Hz0)
jiaodu_0=angle(Hz0)*180/pi


t=0;
for tt=-6:0

ttt=0;
ar=[-3.82704 -3.56608 -3.30512 -3.04416 -2.7832 -2.52224 -2.26128 -2.00032 -1.73936 -1.4784 -1.21744 -0.95648 -0.69552 -0.43456 -0.1736 0.08736 0.34832 0.60928 0.87024 1.1312 1.39216];

w=power(10,(ar-tt));
for j=1:21

for ii=1:3
    bb(ii)=sqrt(i*w(j)*muu*dert(ii)); %第ii層波數(shù)
    dd=sqrt(x.^2+bb(ii).^2); %Uj
    uu(ii)=dd;
end

mmm=coth(uu(1)*h(1)+acoth(uu(1)/uu(2)*coth(uu(2)*h(2)+acoth(uu(2)/uu(3)*1))));

%f=x.*(1+exp(-(jie+fa).*x).*(x.*mmm-uu(1))./(x.*mmm+uu(1)));
f=-i*w(j)*x.*exp(x.*(-jie_z)).*(1+(x.*mmm-uu(1))./(x.*mmm+uu(1)));%電性源半航空


sum=0;
for n=1:47
   lannn=a+(n-1).*s;
   lannni=(1./r).*power(10,lannn); %被積變量
   ff=subs(f,x,lannni);
   z=ff.*han(n); %被積函數(shù)與漢克爾函數(shù)系數(shù)相乘
   sum=sum+z;
end

g=sum./r.^2;


rr_1=sqrt((jie_x-xx)^2+(jie_y).^2); % A邊 

gg_1=subs(g,r,rr_1); % A邊的積分函數(shù)






%AA的simpsion積分
a_1=-bianchang_1;%input('a_1='); %邊長(zhǎng)
b_1=bianchang_1;%input('b_1='); %邊長(zhǎng)
snn2_1=IntGauss('gg_1',a_1,b_1,7)
AA=(jie_y).*snn2_1 


%I=input('I=');
%I=-1;

ttt=ttt+1
Hz(ttt)=I.*(AA)./(4*pi)
Hz_1(ttt)=abs(Hz(ttt))

jiaodu(ttt)=angle(Hz(ttt))*180/pi
end
Hz
t=t+1
fr=[0.000349998 -0.000418371 0.000772828 -0.000171356 0.001022172 0.000897638 0.002208974 0.003844944 0.00680904 0.013029162 0.022661391 0.042972904 0.075423603 0.139346367 0.234486236 0.366178323 0.284615486 -0.235691746 0.046994188 -0.005901946 0.000570165];
ffff(t)=real(Hz)*fr'
aaa(t)=tt
end
ffff
Hz0
ff1=Hz0-ffff
t1=power(10,aaa);
loglog(t1,ff1)
xlabel('t/（s）')
ylabel('Hz/（A/m）')