Numerical-Methods/PlottingBonddistance.c at main · foxmaestro753/Numerical-Methods · GitHub

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
#include <stdio.h>
#include <math.h>

#define EP 8.854187817e-12
#define PI 3.14159
#define Q 1.602e-19
#define V 1.74637253e-16
#define R 0.33e-10

double poten(double x)
{
    double u;
    u = ((Q * Q) / (4 * PI * EP * x)) - (V * exp((-x) / R));
    return u;
}

double potend(double x)
{
    double g;
    double h = 1e-13;
    g = (poten(x + h) - poten(x)) / h;
    return g;
}

int main()
{

    FILE *ptr;

    ptr = fopen("bond1.txt", "w");

    for (double i = 0.5e-12; i <= 1e-10;)
    {
        fprintf(ptr, "%0.18lf %0.50lf\n", i, poten(i));
        i = i + 1e-14;
    }

    ptr = fopen("bond2.txt", "w");

    for (double i = 0.5e-12; i <= 1e-10;)
    {
        fprintf(ptr, "%0.18lf %0.50lf\n", i, potend(i));
        i = i + 1e-14;
    }

    double a, b, d, e;

    a = 1e-14;
    b = 2e-11;

    while (1)
    {

        d = (a + b) / 2;

        if (potend(a) * potend(d) < 0)
        {
            b = d;
        }

        if (potend(b) * potend(d) < 0)
        {
            a = d;
        }

        e = (a + b) / 2;

        if (potend(a) * potend(e) < 0)
        {
            b = e;
        }

        if (potend(b) * potend(e) < 0)
        {
            a = e;
        }

        if (fabs(d-e) < 1e-12)
        {
            break;
        }
    }
    printf("%e", e);
}