/*! \file  MeshProfileIntersectionx.c
 */

#include "./MeshProfileIntersectionx.h"

void MeshProfileIntersectionx( double** psegments, int* pnumsegs, int* index, double* x, double* y, int nel, int nods,  Contour<IssmPDouble>** contours,int numcontours){

	int i,j,k;

	/*Contour:*/
	Contour<IssmPDouble>* contouri=NULL;
	int      numnodes;
	double*  xc=NULL;
	double*  yc=NULL;

	/*output: */
	double* segments=NULL;
	int     numsegs;

	/*intermediary: */
	double** allsegments=NULL;
	double* segmentsi=NULL;
	int*    allnumsegs=NULL;
	int     numsegsi;
	int     count;

	/*Allocate: */
	allsegments=xNew<double*>(numcontours);
	allnumsegs=xNew<int>(numcontours);

	/*Loop through all contours: */
	for (i=0;i<numcontours;i++){

		/*retrieve contour info: */
		contouri=*(contours+i);
		numnodes=contouri->nods;
		xc=contouri->x;
		yc=contouri->y;

		/*determine segmentsi and numsegsi for this contour and the mesh intersection: */
		MeshSegmentsIntersection(&segmentsi,&numsegsi,index,x,y,nel,nods,xc,yc,numnodes);

		/*save segmentsi: */
		allsegments[i]=segmentsi;
		allnumsegs[i]=numsegsi;
	}

	/*total number of segments? */
	numsegs=0;
	for(i=0;i<numcontours;i++)numsegs+=allnumsegs[i];

	/*Out of all segments, create one common array of segments: */
	segments=xNew<double>(5*numsegs);
	count=0;
	for(i=0;i<numcontours;i++){

		segmentsi=allsegments[i];
		numsegsi=allnumsegs[i];

		for(j=0;j<numsegsi;j++){
			for(k=0;k<5;k++){
				*(segments+count*5+k)=*(segmentsi+j*5+k);
			}
			count++;
		}
	}

	/*Assign output pointers:*/
	*psegments=segments;
	*pnumsegs=numsegs;
}
