%單位點源對應全空間格林函數（時間域）
clc
clear
close?all
u=4*pi*10^（-7）;%取為真空磁導率
o=0.01;?%電導率S/m
%r確定t變化
r=100;%觀測點距點源距，以電源為坐標原點
v=-7:0.05:2;
t=10.^v;
Gt=（（u*o）^.5）./（8*（pi*t）.^1.5）.*exp（-u*o*r^2./（4.*t））;
set（gcf‘Position‘[160?100?1000?500]）;
set（gca‘Position‘[.13?.17?.8?.74]）;
?subplot（121）
loglog（t*1000Gt‘-*r‘‘MarkerSize‘3）;
title（‘g-t‘‘fontsize‘13‘Fontname‘‘Times?New?Roman‘）;
xlabel（‘t（ms）‘‘fontsize‘13）;
ylabel（‘g‘‘fontsize‘13）
axis（[10^-4?10^1?10^-4?10^1]）;

%t確定r變化
t=1;?%ms
x=0;
y=0;
z=1:3:100;
r=（x.^2+y.^2+z.^2）.^（0.5）;%觀測點距點源距，以電源為坐標原點
Gr=（（u*o）^.5）*exp（-u*o*r.^2/（4*t））/（8*（pi*t）^（1.5））*（t>=0）;
subplot（122）
plot（rGr‘-*r‘）;
title（‘g-r‘‘fontsize‘13‘Fontname‘‘Times?New?Roman‘）
xlabel（‘r（m）‘‘fontsize‘13）
ylabel（‘g‘‘fontsize‘13）
grid;

%三維輸出
u=4*pi*10^（-7）;%取為真空磁導率
o=0.01;?%電導率S/m
r=1:0.5:91;%觀測點距點源距，以電源為坐標原點
v=-7:0.05:2;
t=10.^v;
[rt]=meshgrid（rt）;
G=（（u*o）^.5）./（8*（pi*t）.^1.5）.*exp（-u*o*r^2./（4.*t））;
figure（‘Numbertitle‘‘off‘‘name‘‘3D-g（tr）‘）;