Здравствуйте, Посетитель - 383763!
Идея алгоритма такова:
Мы не будем искать точных точек на кривой, т.к. их там нет. Не будем строить и сплайны, нет необходимости.
Учтем, что на данном отрезке кривая практически вертикальная, т.е. значения X у соседних точек практически одинаковые.
Мы будем искать точки,
ближайшие к середине!
Итак, делаем следующее:
1) Ищем среднее значение по Y, для этого ищем минимальное и максимальное.
2) Ищем первую точку, для которой расстояние до теоретической средней минимальное.
Ищем до тех пор, пока идет уменьшение.
Чтобы исключить возможную "неправильную" кривизну на краях кривой,
не будем рассматривать точки, для которых длина больше некоторого значения, скажем 100.
3) Перед тем, как найти вторую точку, найдем самую удаленную точку.
4) Повторим этап 2), чтобы найти вторую точку.
Я с Fortran-ом не работаю. Вот Вам реализация на C++
Конечно, программа далека от совершенства, призвана только показать идею алгоритма.
Для упрощения, я убрал из файла первую строку с названиями полей, чтобы были только данные.
Программа дала результат 830.259
CalcF.cpp
[code h=207]#include <iostream>
#include <malloc.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
using namespace std;
typedef struct POINT
{
double x;
double y;
}POINT;
//Загрузка всех точек в массив
//возвращает заполненный массив и количество точек
int ReadPoints(char* fName, POINT** points)
{
char line[256];
char *pline;
POINT one;
int max = 1000;
int count;
*points = (POINT*)malloc((sizeof POINT)*max);
FILE * pFile = fopen( fName, "r" );
if ( !pFile )
return 0;
for (count=0 ;0!=fgets( line, 256, pFile );)
{
one.x = atof(line);
pline = strchr(line, ' ');
one.y = atof(pline);
(*points)[count++] = one;
if (count >= max)
{
max += 1000;
*points = (POINT*)realloc(*points, (sizeof POINT)*max);
}
}
fclose( pFile );
return count;
}
//расчет коэрцитивной силы
//параметры: число точек и массив с точками
double Calc(int count, POINT *points)
{
double ymin, ymax, ymid, ycurr;
double x1, x2;
int i;
//найдем максимальный и минимальный Y
ymax = -10000000; //max
ymin = 10000000; //min
for (i=0; i<count; i++)
{
if (points[i].y > ymax)
ymax = points[i].y;
if (points[i].y < ymin)
ymin = points[i].y;
}
//теоритическое среднее
ymid = (ymax - ymin)/2;
//ищем первую точку
//ycurr - текущее, ближайшее к ymid Y
for (ycurr=100000,i=0; i<count; i++)
{
//идем по точкам, пока идет уменьшение
if (fabs(ymid - points[i].y) < fabs(ymid - ycurr))
{
//запоминаем
ycurr = points[i].y;
x1 = points[i].x;
} //учтем возможную "неправильность" данных на краях
else if (fabs(ymid - points[i].y) < 100)
//выходим из поиска, если получаем вблизи точки увеличение расстояния
break;
}
//"разбегаемся" :), находим наиболее удаленную точку
for (; i<count; i++)
{
if (fabs(ymid - points[i].y) > fabs(ymid - ycurr))
ycurr = points[i].y;
else
break;
}
//и "с разбега", аналогично предыдущему, находим вторую точку
for (; i<count; i++)
{
if (fabs(ymid - points[i].y) < fabs(ymid - ycurr))
{
ycurr = points[i].y;
x2 = points[i].x;
}
else if (fabs(ymid - points[i].y) < 100)
break;
}
//осталось только посчитать модуль разности
return (fabs(x2 - x1));
}
int main(int argc, char *argv[ ])
{
system("chcp 1251 >> nul");
int count;
POINT * points = 0;
if (argc > 1)
{
if (0 != (count = ReadPoints(argv[1], &points)))
{
cout << "Коэрцитивная сила = " << Calc(count, points) << endl;
free(points);
}
else
cout << "Ошибка чтения файла" << endl;
}
else
cout << "Введите: CalcF "<имя файла>" << endl;
system("pause");
return 0;
}
[/code]