## Code for Figure 5.2

### main.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% Copyright (C) 2001, James B. Rawlings and John G. Ekerdt
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License as
% published by the Free Software Foundation; either version 2, or (at
% your option) any later version.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; see the file COPYING.  If not, write to
% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
% MA 02111-1307, USA.

initdata ();

an1=0;
en1=3;

an2=0;
en2=3;

n1 = 31;
n2 = 31;

step1 = (en1 - an1) / (n1 - 1);
step2 = (en2 - an2) / (n2 - 1);

d3table = zeros (n2, 3*n1);
table = zeros (n1, n2);

idx = [1,2,3];
for p1 = 1:n1
r1 = (p1-1)*step1 + an1;
for p2 = 1:n2
r2 = (p2-1)*step2 + an2;
tmp = Vau (r1, r2);
d3table(p2,idx) = [r1, r2, tmp];
table(p1,p2) = tmp;
end
idx += 3;
end

[s, V, d3trace, num] = inittrace (en1, en2);

[fid, errmsg] = fopen ('hfs.dat', 'w');
fprintf(fid, '% [min]\n');
for i = 1:num
j = i*3-2;
fprintf (fid, '%f %f %f\n', d3trace(j:j+2));
end

fprintf(fid, '\n% [surface]\n');
for p2 = 1:n2
fprintf (fid,'%f %f %f\n', d3table(p2,:));
fprintf (fid,'\n');
end
fclose (fid);

```

### Vau.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% Copyright (C) 2001, James B. Rawlings and John G. Ekerdt
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License as
% published by the Free Software Foundation; either version 2, or (at
% your option) any later version.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; see the file COPYING.  If not, write to
% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
% MA 02111-1307, USA.

function retval = Vau(r1,r2)
global D1 D3 B1 B3 r10 r30 R_1 C Si A
r3=r1+r2;
r(1)=r1;
r(2)=r2;
r(3)=r3;
for x = 1:3
E1(x)=D1(x)*(exp(-2*B1(x)*(r(x)-r10(x)))-2*exp(-B1(x)*(r(x)-r10(x))));
E3(x,1)=D3(x)*(exp(-2*B3(x)*(r(x)-r30(x)))+2*exp(-B3(x)*(r(x)-r30(x))));
E3(x,2)=C(x)*(r(x)+A(x))*exp(-Si(x)*r(x));
if(r(x) <= R_1(x))
y=1;
else
y=2;
end
Q(x)=(E1(x)+E3(x,y))/2;
J(x)=(E1(x)-E3(x,y))/2;
end
retval=Q(1)+Q(2)+Q(3)-sqrt(0.5*((J(1)-J(2))*(J(1)-J(2))+(J(2)-J(3))*(J(2)-J(3))+(J(1)-J(3))*(J(1)-J(3))));
end%function

```

### initdata.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% Copyright (C) 2001, James B. Rawlings and John G. Ekerdt
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License as
% published by the Free Software Foundation; either version 2, or (at
% your option) any later version.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; see the file COPYING.  If not, write to
% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
% MA 02111-1307, USA.

function initdata ()

global D1 D3 B1 B3 r10 r30 R_1 C Si A

D1 = [129.3; 94.8; 129.3];
B1 = [2.316; 2.047; 2.316];
r10 = [0.931; 0.753; 0.931];
D3 = [51.83; 34.56; 51.83];
B3 = [2.229; 1.522; 2.229];
r30 = [0.931; 0.753; 0.931];
R_1 = [1.217; 0.443; 1.217];
C = [1174.9; 388.8; 1174.9];
Si = [3.090; 1.929; 3.090];
A = [1.315; 0.762; 1.315];

end%function

```

### inittrace.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% Copyright (C) 2001, James B. Rawlings and John G. Ekerdt
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License as
% published by the Free Software Foundation; either version 2, or (at
% your option) any later version.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; see the file COPYING.  If not, write to
% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
% MA 02111-1307, USA.

function [s, V, d3trace, num] = inittrace (en1, en2)

global D1 D3 B1 B3 r10 r30 R_1 C Si A

R1sar = 1.29781;
R2sar = 0.80363;
numr = 30;

R1r = linspace (R1sar, en1, numr)';

opts = odeset ('AbsTol', sqrt (eps), 'RelTol', sqrt (eps));
[tsolver, R2r] = ode15s (@rder, R1r, R2sar, opts);

R1sal = 1.297;
R2sal = 0.80363;
numl = 1000;

R2l = linspace (R2sal, en2, numl)';

opts  =  odeset ('AbsTol', sqrt (eps), 'RelTol', sqrt (eps));
[tsolver, R1l] = ode15s (@lder, R2l, R1sal, opts);

R1rn = zeros (numr, 1);
R2rn = zeros (numr, 1);
num = numl+numr;

for i = 1:numr
R1rn(i) = R1r(numr+1-i);
R2rn(i) = R2r(numr+1-i);
end

R1 = [R1rn; R1l];
R2 = [R2rn; R2l];
trtable = [R1, R2];

s = zeros (num, 1);
V = zeros (num, 1);
D = zeros (num, 1);

s(1) = 0;
V(1) = Vau (R1(1), R2(1));
D(1) = 1;

for i = 2:num
D(i) = i;
s(i) = sqrt ((R1(i)-R1(i-1))^2 + (R2(i)-R2(i-1))^2) + s(i-1);
V(i) = Vau (R1(i), R2(i));
end

pot = [s, V, R1, R2, D];

for i = 1:num
j = i*3-2;
d3trace(j) = R1(i);
d3trace(j+1) = R2(i);
d3trace(j+2) = V(i);
end

d3trace = d3trace';

% this seems to be unused.
% for a = 1:3
%   for i = 1:num
%    j = i*3-2;
%    d3sptrace(a,j) = R1(i)+(a-2)*0.1;
%    d3sptrace(a,j+1) = R2(i);
%    d3sptrace(a,j+2) = -100+(a-2);
%   end
% end

end%function

```

### lder.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% Copyright (C) 2001, James B. Rawlings and John G. Ekerdt
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License as
% published by the Free Software Foundation; either version 2, or (at
% your option) any later version.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; see the file COPYING.  If not, write to
% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
% MA 02111-1307, USA.

function retval = lder (r2, r1)

global D1 D3 B1 B3 r10 r30 R_1 C Si A

r3 = r1+r2;
r(1) = r1;
r(2) = r2;
r(3) = r3;

for x = 1:3
E1(x) = D1(x)*(exp(-2*B1(x)*(r(x)-r10(x)))-2*exp(-B1(x)*(r(x)-r10(x))));
E3(x,1) = D3(x)*(exp(-2*B3(x)*(r(x)-r30(x)))+2*exp(-B3(x)*(r(x)-r30(x))));
E3(x,2) = C(x)*(r(x)+A(x))*exp(-Si(x)*r(x));
if(r(x) <= R_1(x))
y = 1;
else
y = 2;
end
Q(x) = (E1(x)+E3(x,y))/2;
J(x) = (E1(x)-E3(x,y))/2;
end

%dQ(a,b) means dQ(a)/dr(b)
for i = 1:3
dE1(i,i) = -2*B1(i)*D1(i)*(exp(-2*B1(i)*(r(i)-r10(i)))-exp(-B1(i)*(r(i)-r10(i))));
if(r(i) <= R_1(i))
dE3(i,i) = -2*B3(i)*D3(i)*(exp(-2*B3(i)*(r(i)-r30(i)))+exp(-B3(i)*(r(i)-r30(i))));
else
dE3(i,i) = C(i)*exp(-Si(i)*r(i))*(1-Si(i)*(r(i)+A(i)));
end
end
dE1(1,2) = 0;
dE3(1,2) = 0;
dE1(2,1) = 0;
dE3(2,1) = 0;
dE1(3,1) = dE1(3,3);
dE3(3,1) = dE3(3,3);
dE1(3,2) = dE1(3,3);
dE3(3,2) = dE3(3,3);
dE1(1,3) = 1000;
dE3(1,3) = 1000;
dE1(2,3) = 1000;
dE3(2,3) = 1000;
dE1(3,3) = 1000;
dE3(3,3) = 1000;
for i = 1:3
for j = 1:3
dQ(i,j) = (dE1(i,j)+dE3(i,j))/2;
dJ(i,j) = (dE1(i,j)-dE3(i,j))/2;
end
end

dV(1) = dQ(1,1)+dQ(3,1)-0.5/sqrt(0.5*((J(1)-J(2))*(J(1)-J(2))+(J(2)-J(3))*(J(2)-J(3))+(J(1)-J(3))*(J(1)-J(3))))*(J(1)*(2*dJ(1,1)-dJ(3,1))+J(2)*(-dJ(1,1)-dJ(3,1))+J(3)*(2*dJ(3,1)-dJ(1,1)));

dV(2) = dQ(2,2)+dQ(3,2)-0.5/sqrt(0.5*((J(1)-J(2))*(J(1)-J(2))+(J(2)-J(3))*(J(2)-J(3))+(J(1)-J(3))*(J(1)-J(3))))*(J(1)*(-dJ(2,2)-dJ(3,2))+J(2)*(2*dJ(2,2)-dJ(3,2))+J(3)*(2*dJ(3,2)-dJ(2,2)));

retval = dV(1)/dV(2);
end%function

```

### rder.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% Copyright (C) 2001, James B. Rawlings and John G. Ekerdt
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License as
% published by the Free Software Foundation; either version 2, or (at
% your option) any later version.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; see the file COPYING.  If not, write to
% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
% MA 02111-1307, USA.

function retval = rder (r1, r2)

global D1 D3 B1 B3 r10 r30 R_1 C Si A

retval = 1 / lder (r2, r1, D1, D3, B1, B3, r10, r30, R_1, C, Si, A);

end%function

```