diff --git a/src/init.c b/src/init.c index 038ebca..1e42587 100644 --- a/src/init.c +++ b/src/init.c @@ -100,8 +100,45 @@ double chi(double b){ return find_root_5ord(y_previous[0], y[0], y_previous[1], y[1], phi_previous, phi); } -static void refine_interval(double xl, double xr, double yl, double yr, SampleData *sample){} +void sample_push(SampleData *sample, double x, double y){ + if (sample->size + 1 >= sample->capacity) { + sample->x = realloc(sample->x, sizeof(double) * (sample->capacity)*2); + sample->y = realloc(sample->y, sizeof(double) * (sample->capacity)*2); + } + sample->x[sample->size] = x; + sample->y[sample->size] = y; + sample->size += 1; +} + +static void refine_interval(double xl, double xr, double yl, double yr, double rela_err_limit, SampleData *sample){ + double xm = (xl + yl)/2; + double ym = chi(xm); + sample_push(sample, xm, ym); + if (fabs((ym-(yl+yr)/2)/ym) >= rela_err_limit) { + refine_interval(xl, xm, yl, ym, rela_err_limit, sample); + refine_interval(xm, xr, ym, yr, rela_err_limit, sample); + } + return ; +} static int init(double bmax, double rela_err_limit){ + SampleData sample; + sample.capacity = 2; + sample.size=2; + sample.x = malloc(sizeof(double) * 2); + sample.y = malloc(sizeof(double) * 2); + sample.x[0] = bmin; + sample.x[1] = bmax; + sample.y[0] = chi(bmin); + sample.y[1] = chi(bmax); + + refine_interval(bmin, bmax, sample.y[0], sample.y[1], rela_err_limit, &sample); + + //sort + + //output + + free(sample.x); + free(sample.y); return 0; }