#include <stdlib.h>
#include <unistd.h>
#include <stdio.h>
#include <string.h>
#include <signal.h>
#include <pthread.h>
#include <sched.h>
#include <math.h>
#include <sys/time.h>
#include <fcntl.h>
#include <unistd.h>

typedef struct
{
	pthread_t tid;
	int first;
	int last;
} ThreadRecord;

ThreadRecord *threads;

double *vec_next ; //массив элементов в следующий момент времени
double *vec ; //массив элементов в текущий момент времени
double *p ;

pthread_barrier_t barr1, barr2;

struct timeval tv1,tv2,dtv;
struct timezone tz;

double PI = 3.1415926535; //пи
double R = 1000; //сопротивление
double C = 1e-4; //ёмкость
double IA = 2 ; //амплитуда силы тока
int done = 0; //признак выполнения
int time_moment = 0; //момент времени
int num_of_elements = 100; //количество элементов(по умолчанию)
int num_of_threads = 8; //количество потоков(по умолчанию)
int time_intervals = 20 ; //количество временных интервалов(по умолчанию)
double step ; //временной шаг

void print_nodes(double *vec_next, FILE **gp)
{ //визуализация работы
	int j ;
	//fprintf (*gp, "plot [0:%d][0:%lf]'-' with line title 'sin(x %+#.3g)'\n",num_of_elements,IA*R + 10 , -j*0.1);
	for ( j = 0 ; j < num_of_elements; j ++)
	{
		fprintf (*gp, " %# -15g\n", vec_next[ j]);
	}
	fprintf (*gp, "\n\n"); fflush (*gp); usleep (10000);
}

void * node_method(void *arg_p)
{ //вычисление в узлах явным методом Эйлера
	int i ;
	ThreadRecord *thr;
	thr = (ThreadRecord *)arg_p;
	while(!done)
	{
		pthread_barrier_wait(&barr1);
		for ( i = thr->first; i < thr->last; i ++)
		{
			if(i && (i < num_of_elements - 1 ) &&(time_moment + 1 < (int)((double)time_intervals/step) ))
			{
				vec_next[ i] = vec[i] + step*(vec[ i - 1] - 2 * vec[i] + vec[i + 1])/(R*C);
			}
		}
	pthread_barrier_wait(&barr2);
	}
}

int main(int argc, char *argv[], char *env[])
{
	step = R*C/10;
	FILE *gp = fopen("data_test.txt", "w"); //файл для визуализации
	int i, j;
	if(argc >= 4)
	{ //получение данных из командной строки
		num_of_elements = atoi(argv[1]);
		num_of_threads = atoi(argv[2]);
		time_intervals = atoi(argv[3]);
		if(num_of_elements < num_of_threads)
			num_of_threads = num_of_elements;
	}
	vec_next = (double *)malloc( num_of_elements * sizeof(double ));
	vec = (double *)malloc(num_of_elements * sizeof(double ));
	for ( i = 0; i < num_of_elements; i ++)
	{
		vec[i] = 0 ;
	}
	pthread_attr_t pattr;
	pthread_attr_init (&pattr);
	pthread_attr_setscope (&pattr, PTHREAD_SCOPE_SYSTEM);
	pthread_attr_setdetachstate (&pattr,PTHREAD_CREATE_JOINABLE);
	threads = (ThreadRecord *)malloc(num_of_threads * sizeof(ThreadRecord));
	pthread_barrier_init(&barr1, NULL, num_of_threads + 1);
	pthread_barrier_init(&barr2, NULL, num_of_threads + 1);
	int d = (num_of_elements%num_of_threads) ? 1 : 0 ;
	int len = num_of_elements/num_of_threads + d ;
	gettimeofday(&tv1, &tz);
	for ( i = 0; i < num_of_threads; i ++)
	{
		threads[i].first = i * len ;
		threads[i].last = ( i + 1 != num_of_threads) ? i * len + len : num_of_elements ;
		if(pthread_create(&(threads[i].tid),&pattr, node_method, (void *)&(threads[i])))
			perror("pthread create");
	}
	vec_next[0] = 0;
	vec_next[num_of_elements - 1] = 0;
	pthread_barrier_wait(&barr1);
	for ( i = 1 ; i <= (int)((double)time_intervals/step); i ++)
	{ //явный метод Эйлера. начало
		pthread_barrier_wait(&barr2);
		time_moment = i ;
		if (time_moment + 1 < (int)((double)time_intervals/(step * 100)) )
		{ //краевые значения для одной десятой времени работы
			vec_next[0 ] = R*IA ;
			vec_next[num_of_elements - 1] = R*IA ;
		}
		else if(time_moment + 1 < (int)((double)time_intervals/(step )))
		{ //краевые значения для остального времени
			vec_next[ 0 ] = vec[1];
			vec_next[ num_of_elements - 1] = vec[ num_of_elements - 2];
		}
		if (argc >= 5)
		{
			print_nodes(vec_next, &gp) ; //вывод результатов визуализации
		}
		p = vec ;
		vec = vec_next ;
		vec_next = p ;
		pthread_barrier_wait(&barr1);
	}done = 1 ;
	gettimeofday(&tv2, &tz);
	dtv.tv_sec= tv2.tv_sec -tv1.tv_sec;
	dtv.tv_usec=tv2.tv_usec-tv1.tv_usec;
	if(dtv.tv_usec<0)
	{
		dtv.tv_sec--; dtv.tv_usec+=1000000;
	}
	printf("Time:%ld.%ld\n",(dtv.tv_sec*1000000+dtv.tv_usec)/1000000,(dtv.tv_sec*1000000+dtv.tv_usec)%1000000); //время работы
}