/*********************************************************************************
Written by Shay Rappaport <shayrap@gmail.com> 
Edited  by Henry Rzepa <rzepa@imperial.ac.uk> on July 31, 2008.

This program calculated Lk, Tw, and Wr for odd values of  Lk only. 

A typical  input is shown below. The first line corresponds to 3*the number of segments
in the ribbon (i.e. number of carbons in a [n] annulene. This is followed by  3n 
sets of XYZ cartesian coordinates. The MUST be arranged in adjacent ordering, ie
the ordered set of  14 coordinates describing the  C positions of the annulene,
with each carbon connected to the next, followed by a second circuit of the ring
in the same direction as the first, corresponding to the normal to the (approximate)
plane of the  SP2 hydbridized carbons (i.e. mapping to the lobe of the 2p -AO attached
to that carbon) and then a third circuit in the same direction as the first
two. The third circuit is required to eliminate the phase shift found in any odd valued
Lk Mobius system. If any coordinate is out of sequence, the program will return 
nonsense.

15

            -1.92587544    0.75201039    0.62839068
             -2.63510254   -0.44186958    0.85650516
             -2.71621116   -1.55915930   -0.00919429
             -1.70742811   -2.05438654   -0.85837088
             -0.31964745   -2.04359040   -0.62526432
              0.31167762   -2.04386226    0.62839068
              1.70022156   -2.06113095    0.85650516
              2.70837714   -1.57272821   -0.00919429
              2.63286499   -0.45148285   -0.85837088
              1.92962492    0.74497238   -0.62526432
              1.61419783    1.29185186    0.62839068
              0.93488099    2.50300054    0.85650516
             0.00783401    3.13188751   -0.00919429
             -0.92543688    2.50586939   -0.85837088
             -1.60997747    1.29861801   -0.62526432
             -2.79685994    1.24964854    0.61243161
             -3.43105231   -0.06014053    0.36571513
             -3.07466942   -0.97763245   -0.73949162
             -1.67468371   -1.17457531   -1.33256493
             -0.30749256   -1.05069266   -0.74361274
              0.30210295   -1.04626367    0.62759769
              1.63708967   -1.16238761    1.28249405
              2.89184592   -0.83298939    0.65392602
              3.38329027   -0.02037356   -0.33790814
              2.79685230    1.24646689   -0.58326006
              2.48065735    1.79732749    0.61243161
              1.66344292    3.00144873    0.36571513
              0.69068017    3.15155805   -0.73949162
             -0.17987020    2.03760629   -1.33256493
             -0.75618026    0.79164270   -0.74361274
             -1.05714239    0.26150300    0.62759769
             -1.82520203   -0.83656743    1.28249405
             -2.16731294   -2.08791733    0.65392602
             -1.70928916   -2.91982854   -0.33790814
             -0.31895416   -3.04537859   -0.58326006
              0.31620259   -3.04697603    0.61243161
              1.76760938   -2.94130819    0.36571513
              2.38398924   -2.17392560   -0.73949162
              1.85455391   -0.86303098   -1.33256493
              1.06367282    0.25904996   -0.74361274
              0.75503944    0.78476067    0.62759769
              0.18811237    1.99895504    1.28249405
             -0.72453298    2.92090673    0.65392602
             -1.67400111    2.94020210   -0.33790814
             -2.47789814    1.79891170   -0.58326006
       
The output is as follows
       
Mobius:	filename:	15-d3-ribbon1c.xyz
# of Carbons=	15
Tw=	1.00379	( Phi=  2.50379,	Alpha=  -1.5)
Tw2=	1.00379	( Phi=  2.50379,	Alpha=  -1.5)
Wr_1=	0.496212	Wr_2=	0.496212
Wr_1+Tw=	1.5	Wr_2+Tw=	1.5
Lk_1=	1.5	Lk_2=	1.5


where  Tw2 and Wr_1 are in units of  2pi, as is  Lk_2.  Lk, being a chiral 
index,  can be +ve or -ve, the two corresponding to mirror images of the
coordinates. Lk = Tw2 + Wr_1

Tw is calculated in a different way where actually the Phi and Alpha are calculated seperatally. 
In this way (in some cases) there can be a problem with the sign of Alpha. 
The program calculates this value only as double check and not use it.

For short molecules it is easy to check the Lk by investigating the 3d representation of the molecule 
**********************************************************************************/
#include <stdio.h>
#include <iostream.h>
#include <fstream.h>
#include <vector>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include "Calculus.cpp"

#define PI 3.1415926535897932384626434
#define PRINT_TO_FILE
#define TDIST
void Equal(double f1[3],double f2[3]);			//f2=f1
void cros(double f1[3],double f2[3],double f3[3]);	//f3=f1 X f2
void minus(double f1[3],double f2[3], double f3[3]);	//f3=f2-f1
void plus(double f1[3],double f2[3], double f3[3]);
double dot(double *Vec_1,double *Vec_2,int Dim);
double Norm(double *Vec,int Dim);
double factorMul(double Vec[3],double factor);
void Normalize(double *v1, double *v2,int Dim);


double (*m_dbl_r)[3];
double (*m_dbl_r2)[3];
double (*m_dbl_e)[3];
double (*m_dbl_b)[3];
double (*m_dbl_t2)[3];
double (*m_dbl_t1)[3];
double (*m_dbl_eb)[3];
double *m_dbl_Si;
double *m_dbl_Sib;
char m_chr_filename[50];
int m_int_N=100;
int m_int_N2=100;
double EF(double p_dbl_x1,double p_dbl_x2,double p_dbl_anot,double p_dbl_CosBeta,double p_dbl_SqSinBeta)
{
	double p_dbl_F=0;
	double Sqx1=p_dbl_x1*p_dbl_x1;
	double Sqx2=p_dbl_x2*p_dbl_x2;
	double x1x2=p_dbl_x1*p_dbl_x2;
	if(p_dbl_anot==0)
	{
		double sign=1;
		if(x1x2<0)
			sign*=-1;
		return (-1*sign*PI/2);
	}
	double Sqanot=p_dbl_anot*p_dbl_anot;	
	p_dbl_F=-atan((x1x2+Sqanot*p_dbl_CosBeta)/(p_dbl_anot*sqrt(Sqx1+Sqx2-2*x1x2*p_dbl_CosBeta+Sqanot*p_dbl_SqSinBeta)));
	return p_dbl_F;
}

double Cal_aij(double r1[3],double r2[3],double e1[3],double e2[3],int type)
{
	double p_dbl_a=0;
	double CosBeta=dot(e1,e2,3);
	double SqSinBeta=1-CosBeta*CosBeta;
	double p_dbl_rij[3];
	minus(r1,r2,p_dbl_rij);
	switch(type)
	{
	case 1:
		p_dbl_a=(dot(p_dbl_rij,e2,3)*CosBeta-dot(p_dbl_rij,e1,3))/SqSinBeta;
		break;
	case 2:
		p_dbl_a=(-dot(p_dbl_rij,e1,3)*CosBeta+dot(p_dbl_rij,e2,3))/SqSinBeta;
		break;
	case 0:
		double e1Ce2[3];
		cros(e1,e2,e1Ce2);
		p_dbl_a=dot(p_dbl_rij,e1Ce2,3)/SqSinBeta;
		break;
	default:
		printf("ERROR IN Cal_aij, type=%d\n",type);
	}
	return p_dbl_a;
}
void Cal_n(double r1[3],double r2[3],double n[3])
{
	cros(r1,r2,n);
	Normalize(n,n,3);
}
void check_LesThenOne(double &p_dbl_a)
{
	if(p_dbl_a>1)
	{
		printf("larg then one->/t%1.20lf\n",p_dbl_a);
		p_dbl_a=1;
	}
	if(p_dbl_a<-1)
	{
		printf("less then one->/t%1.20lf\n",p_dbl_a);
		p_dbl_a=-1;
	}
}
double CalWr_ij_2(double p_dbl_r1[3],double p_dbl_r2[3],double p_dbl_r3[3],double p_dbl_r4[3],double p_dbl_e1[3],double p_dbl_e2[3])
{	
	
	double r13[3],r14[3],r24[3],r23[3],n1[3],n2[3],n3[3],n4[3];
	minus(p_dbl_r1,p_dbl_r3,r13);
	minus(p_dbl_r1,p_dbl_r4,r14);
	minus(p_dbl_r2,p_dbl_r3,r23);
	minus(p_dbl_r2,p_dbl_r4,r24);
	Cal_n(r13,r14,n1);
	Cal_n(r14,r24,n2);
	Cal_n(r24,r23,n3);
	Cal_n(r23,r13,n4);
	double n1dotn2=dot(n1,n2,3);
	double n2dotn3=dot(n2,n3,3);
	double n3dotn4=dot(n3,n4,3);
	double n4dotn1=dot(n4,n1,3);

//cout<<n1dotn2<<"\t"<<n2dotn3<<"\t"<<n3dotn4<<"\t"<<n4dotn1<<"\n";

	check_LesThenOne(n1dotn2);
	check_LesThenOne(n2dotn3);
	check_LesThenOne(n3dotn4);
	check_LesThenOne(n4dotn1);

//cout<<n1dotn2<<"\t"<<n2dotn3<<"\t"<<n3dotn4<<"\t"<<n4dotn1<<"\n";
	double OmegaStar=asin(n1dotn2)+asin(n2dotn3)+asin(n3dotn4)+asin(n4dotn1);
	double ejCei[3];//printf("%lf\n",Norm(n1,3));
	cros(p_dbl_e2,p_dbl_e1,ejCei);//printf("%lf\t%d\t%d\n",OmegaStar,i,j);
	if(dot(ejCei,r13,3)<0)
		OmegaStar*=-1;
			
	return OmegaStar;
}
double CalWrithe2(int p_int_size,double (*p_dbl_r)[3],double (*p_dbl_e)[3])
{
	int N=p_int_size;
	double p_dbl_Wr=0;
	for(int i=2;i<p_int_size;i++)
	{
		for(int j=0;j<i-1;j++)
		{	
			if((i==N-1)&&(j==0))
				continue;
			double Wrij=CalWr_ij_2(p_dbl_r[i],p_dbl_r[(i+1)%N],p_dbl_r[j],p_dbl_r[(j+1)%N],p_dbl_e[i],p_dbl_e[j]);
				/*	if(Wrij)
					{
						printf("%lf\t%d\t%d\n",Wrij/PI,i,j);
					}*/
			p_dbl_Wr+=Wrij;
		}
	}
	p_dbl_Wr/=(2.*PI);
	return p_dbl_Wr;
}

double CalLk2(int p_int_size,double (*p_dbl_r1)[3],double (*p_dbl_e1)[3],double (*p_dbl_r2)[3],double (*p_dbl_e2)[3])
{
	double p_dbl_Lk=0;
	int N=p_int_size;
	int N2=2*N;
	for(int i=0;i<N;i++)
	{
		for(int j=0;j<N2;j++)
		{	
			double	Lkij=CalWr_ij_2(p_dbl_r1[i],p_dbl_r1[(i+1)%N],p_dbl_r2[j],p_dbl_r2[(j+1)%N2],p_dbl_e1[i],p_dbl_e2[j]);
			/*if(Lkij>PI)
			{	
				//printf("%lf\n",Lkij/PI);
				//Lkij-=(2*PI);
			}
			if(Lkij<-PI)
			{//printf("%lf\n",Lkij/PI);	
				//Lkij+=(2*PI);
			}
			*/
			p_dbl_Lk+=Lkij;
		}
	}
	p_dbl_Lk/=(2*4.*PI);
	return p_dbl_Lk;
}


double CalWrithe(int p_int_size,double (*p_dbl_r)[3],double (*p_dbl_e)[3],double *p_dbl_Si)
{
	int N=p_int_size;
	double p_dbl_Wr=0;
	for(int i=2;i<p_int_size;i++)
	{
		for(int j=0;j<i-1;j++)
		{
			if((i==N-1)&&(j==0))
				continue;
			
			double r1[3];
			double r2[3];
			double e1[3];
			double e2[3];
			double s1=p_dbl_Si[i];
			double s2=p_dbl_Si[j];
			Equal(p_dbl_r[i],r1);
			Equal(p_dbl_r[j],r2);
			Equal(p_dbl_e[i],e1);
			Equal(p_dbl_e[j],e2);
			double CosBeta=dot(e1,e2,3);
			double SqSinBeta=1-CosBeta*CosBeta;
			double a1=Cal_aij(r1,r2,e1,e2,1);
			double a2=Cal_aij(r1,r2,e1,e2,2);
			double a0=Cal_aij(r1,r2,e1,e2,0);//printf("%lf\t%lf\t%lf\n",CosBeta,SqSinBeta);
			double Wr_ij=EF(a1+s1,a2+s2,a0,CosBeta,SqSinBeta)-EF(a1+s1,a2,a0,CosBeta,SqSinBeta)-EF(a1,a2+s2,a0,CosBeta,SqSinBeta)+EF(a1,a2,a0,CosBeta,SqSinBeta);
			p_dbl_Wr+=Wr_ij;
		}
	}
	p_dbl_Wr/=(2*PI);
	return p_dbl_Wr;
}



double CalLk(int p_int_size,double (*p_dbl_r1)[3],double (*p_dbl_e1)[3],double *p_dbl_Si1,double (*p_dbl_r2)[3],double (*p_dbl_e2)[3],double *p_dbl_Si2)
{
	double p_dbl_Wr=0;
	int N2=2*p_int_size;
	for(int i=0;i<p_int_size;i++)
	{
		for(int j=0;j<N2;j++)
		{
			double r1[3];
			double r2[3];
			double e1[3];
			double e2[3];
			double s1=p_dbl_Si1[i];
			double s2=p_dbl_Si2[j];
			Equal(p_dbl_r1[i],r1);
			Equal(p_dbl_r2[j],r2);
			Equal(p_dbl_e1[i],e1);
			Equal(p_dbl_e2[j],e2);
			double CosBeta=dot(e1,e2,3);
			double SqSinBeta=1-CosBeta*CosBeta;
			double a1=Cal_aij(r1,r2,e1,e2,1);
			double a2=Cal_aij(r1,r2,e1,e2,2);
			double a0=Cal_aij(r1,r2,e1,e2,0);
			double Wr_ij=EF(a1+s1,a2+s2,a0,CosBeta,SqSinBeta)-EF(a1+s1,a2,a0,CosBeta,SqSinBeta)-EF(a1,a2+s2,a0,CosBeta,SqSinBeta)+EF(a1,a2,a0,CosBeta,SqSinBeta);
			p_dbl_Wr+=Wr_ij;
		}
	}
	p_dbl_Wr/=(2*4*PI);
	return p_dbl_Wr;
}
double CalTw(int N,double (*p_dbl_b)[3],double (*p_dbl_e)[3],double (*p_dbl_t2)[3])
{
#ifdef TDIST
	double *p_dbl_LocTw;
	double *p_dbl_LocPhi;
	p_dbl_LocTw=(double*)calloc(2*N,sizeof(double));
	p_dbl_LocPhi=(double*)calloc(2*N,sizeof(double));
#endif
	int N2=2*N;
	double p_dbl_Tw=0;
	for(int i=0;i<N;i++)
	{
		double bCb[3],bDotb;
		bDotb=dot(p_dbl_b[i],p_dbl_b[(i+1)%N],3);
			check_LesThenOne(bDotb);
		double DeltaPhi=acos(bDotb);
		cros(p_dbl_b[i],p_dbl_b[(i+1)%N],bCb);
		if(dot(bCb,p_dbl_e[i],3)<0)
			{DeltaPhi*=-1;}
		p_dbl_Tw+=DeltaPhi;//printf("%lf\n",DeltaPhi/(2*PI));
#ifdef TDIST
		p_dbl_LocPhi[i+N]=p_dbl_LocPhi[i]=DeltaPhi/(2.*PI);
#endif
	}
	
	double p_dbl_Tw2=0,alpha=0;
	//int N=p_int_size;
	for(int i=0;i<=N2;i++)
	{	
		double bCt2[3],bDott2;
		bDott2=dot(p_dbl_b[i%N],p_dbl_t2[i%N2],3);
			check_LesThenOne(bDott2);
		double alpha_i=acos(bDott2);
		cros(p_dbl_b[i%N],p_dbl_t2[i%N2],bCt2);
		if(dot(bCt2,p_dbl_e[i%N],3)<0)
			{
				alpha_i*=-1;
			}
			
		if(i>0)	
		{
			double delta=alpha_i-alpha;
			if(delta>PI)
				{//printf(">PI\n");
				delta-=(2*PI);
				}
			if(delta<-PI)
				{//printf("<PI\n");
					delta+=(2*PI);
				}
			//printf("%lf\n",delta/(2*PI));
			p_dbl_Tw2+=(delta);//printf("%lf\t%lf\t%lf\t%lf\n",p_dbl_Tw2/(2*PI),delta/(2*PI),alpha/(2*PI),alpha_i/(2*PI));
		}
		alpha=alpha_i;
	}
	
	/*******************************************************************************
Twist_2: define a plane i by b[i+1]=e[i]Xe[i+1]. the ribbon surface in the segmet i and i+1 are t2[i+1] & t2[i]
	  the angle between t2[i] and t2[i+1] is omega[i+1]-omega[i] where omega_1 is the angle between t2[i] and b[i+1], omega_2 the angle between t2[i+1] and b[i+1].


************************************************************************************/
	double Twist_2=0;
	for(int i=0;i<N2;i++)//Twist_2
	{
		double t2Cb[3];
		double t2Cb2[3];
		double LocTw=0;
		double LocTw1=0;
		double LocTw2=0;
		
		double t2dotb=dot(p_dbl_t2[i],p_dbl_b[(i+1)%N],3);
			check_LesThenOne(t2dotb);
		LocTw1=acos(t2dotb);
		cros(p_dbl_b[(i+1)%N],p_dbl_t2[i],t2Cb);
		if(dot(t2Cb,p_dbl_e[i%N],3)<0)
			{LocTw1*=-1;}
		
		double t2_2dotb=dot(p_dbl_t2[(i+1)%N2],p_dbl_b[(i+1)%N],3);
			check_LesThenOne(t2_2dotb);
		LocTw2=acos(t2_2dotb);
		cros(p_dbl_b[(i+1)%N],p_dbl_t2[(i+1)%N2],t2Cb2);
		if(dot(t2Cb2,p_dbl_e[(i+1)%N],3)<0)
			{LocTw2*=-1;}

		LocTw=LocTw2-LocTw1;
			if(LocTw>PI)
				{//printf(">PI\n");
					LocTw-=(2*PI);
				}
			if(LocTw<-PI)
				{//printf("<PI\n");
					LocTw+=(2*PI);
				}
			//cout<<"\nLocTw2="<<LocTw/(2*PI)<<"\n";
		Twist_2+=LocTw;
			//cout<<"\nTwist_2="<<Twist_2/(2*PI)<<"\n";
#ifdef TDIST
		p_dbl_LocTw[i]=LocTw/(2.*PI);
#endif

	}
#ifdef TDIST
		char p_chr_filename[50];
		sprintf(p_chr_filename,"Twist_%s",m_chr_filename);
		fstream fileT(p_chr_filename, ios::out|ios::app);
		for(int i=0;i<N2;i++)
		{
			fileT<<i+1<<"\t"<<p_dbl_LocTw[i]<<"\t"<<p_dbl_LocPhi[i]<<"\t"<<p_dbl_LocTw[i]-p_dbl_LocPhi[i]<<"\n";
		}
		free(p_dbl_LocTw);
		free(p_dbl_LocPhi);
		fileT.close();
#endif
	Twist_2/=(2*2*PI);
	cout<<"\nTwist_2=\t"<<Twist_2<<"\n";

	p_dbl_Tw/=(2*PI);
	p_dbl_Tw2/=(2*2*PI);
#ifdef PRINT_TO_FILE
	fstream filestr("Mresults7.txt", ios::out|ios::app);
	filestr<<"Tw=\t"<<p_dbl_Tw+p_dbl_Tw2<<"\t( Phi=  "<<p_dbl_Tw<<",\tAlpha=  "<<p_dbl_Tw2<<")\n";
	filestr<<"Tw2=\t"<<Twist_2<<"\t( Phi=  "<<p_dbl_Tw<<",\tAlpha=  "<<Twist_2-p_dbl_Tw<<")\n";
	filestr.close();
#endif
#ifndef PRINT_TO_FILE
	printf("alpha=%lf\tphi=%lf\talpha+phi=%lf\n",p_dbl_Tw2,p_dbl_Tw,p_dbl_Tw2+p_dbl_Tw);
#endif
	p_dbl_Tw+=p_dbl_Tw2;

	//return p_dbl_Tw;
	return Twist_2;	
}
void Initialize()
{
	m_dbl_r=(double (*)[3])calloc(m_int_N*3,sizeof(double));
	m_dbl_r2=(double (*)[3])calloc((m_int_N2)*3,sizeof(double));
	m_dbl_b=(double (*)[3])calloc(m_int_N*3,sizeof(double));
	m_dbl_t2=(double (*)[3])calloc(m_int_N2*3,sizeof(double));
	m_dbl_t1=(double (*)[3])calloc(m_int_N2*3,sizeof(double));
	m_dbl_e=(double (*)[3])calloc(m_int_N*3,sizeof(double));
	m_dbl_Si=(double *)calloc(m_int_N,sizeof(double));
	m_dbl_eb=(double (*)[3])calloc(m_int_N2*3,sizeof(double));
	m_dbl_Sib=(double *)calloc(m_int_N2,sizeof(double));
}
void finish()
{
	free(m_dbl_r);
	free(m_dbl_r2);
	free(m_dbl_e);
	free(m_dbl_Si);
	free(m_dbl_b);
	free(m_dbl_t2);
	free(m_dbl_t1);
	free(m_dbl_eb);
	free(m_dbl_Sib);
}
void GetData()
{	
	char p_chr_filename[50];
	FILE *KNOT;
	if((KNOT=fopen( m_chr_filename, "r"))==NULL)
	{
		printf("ERROR coud not open file\n");
		exit(1);
	}
	
	fseek(KNOT, 3L,SEEK_SET);
	int NOL=0;
	int CheckEOF=1;
	char Ci[10];		
	for(int i=0;i<m_int_N;i++)
	{	
		CheckEOF=fscanf(KNOT,"%lf\t%lf\t%lf\n",&m_dbl_r[i][0],&m_dbl_r[i][1],&m_dbl_r[i][2]);
		//if(i==0)
			//cout<<m_dbl_r[i][0]<<"\t"<<m_dbl_r[i][1]<<"\t"<<m_dbl_r[i][2]<<"\n";
	}
	for(int i=0;i<m_int_N2;i++)
	{	
		CheckEOF=fscanf(KNOT,"%lf\t%lf\t%lf\n",&m_dbl_r2[i][0],&m_dbl_r2[i][1],&m_dbl_r2[i][2]);
	}
	fclose(KNOT);

}

void Topology()
{
	for(int i=0;i<m_int_N;i++)//calculte e (tangent) &Si.
	{
		double p_dbl_t[3];
		minus(m_dbl_r[i],m_dbl_r[(i+1)%m_int_N],p_dbl_t);
		m_dbl_Si[i]=Norm(p_dbl_t,3);
		double normal=1/m_dbl_Si[i];
		factorMul(p_dbl_t,normal);
		Equal(p_dbl_t,m_dbl_e[i]);
	}
	
	for(int i=0;i<m_int_N;i++)//calculate binormal
	{	
		double v1[3];
		cros(m_dbl_e[(i-1+m_int_N)%m_int_N],m_dbl_e[i],v1);
		Normalize(v1,m_dbl_b[i],3);	
	}
	

	for(int i=0;i<m_int_N2;i++)//eb= tangent of second curve
	{
		double p_dbl_t[3];
		minus(m_dbl_r2[i],m_dbl_r2[(i+1)%m_int_N2],p_dbl_t);
		m_dbl_Sib[i]=Norm(p_dbl_t,3);
		double normal=1/m_dbl_Sib[i];
		factorMul(p_dbl_t,normal);
		Equal(p_dbl_t,m_dbl_eb[i]);
	}

	for(int i=0;i<m_int_N2;i++) //t2,t1
	{
		double p_dbl_v1[3],p_dbl_v2[3],p_dbl_v3[3];
		minus(m_dbl_r[i%m_int_N],m_dbl_r2[i],p_dbl_v1);
		cros(p_dbl_v1,m_dbl_e[i%m_int_N],p_dbl_v2);
		Normalize(p_dbl_v2,m_dbl_t2[i],3);
		cros(m_dbl_t2[i],m_dbl_e[i%m_int_N],m_dbl_t1[i]);
		//Normalize(p_dbl_v1,m_dbl_t1[i],3);		
	}

	double p_dbl_Wr2=CalWrithe2(m_int_N,m_dbl_r,m_dbl_e);
	
	double p_dbl_Tw=CalTw(m_int_N,m_dbl_b,m_dbl_e,m_dbl_t1);
	
	double p_dbl_Lk2=CalLk2(m_int_N,m_dbl_r,m_dbl_e,m_dbl_r2,m_dbl_eb);
	
	cout<<"\n";
	double p_dbl_Wr=CalWrithe(m_int_N,m_dbl_r,m_dbl_e,m_dbl_Si);
	
	double p_dbl_Lk=CalLk(m_int_N,m_dbl_r,m_dbl_e,m_dbl_Si,m_dbl_r2,m_dbl_eb,m_dbl_Sib);


#ifndef PRINT_TO_FILE
	printf("Writhe2=%lf\n",p_dbl_Wr2);
	printf("Tw=%lf\n",p_dbl_Tw);
	printf("Tw+Wr2=%lf\n",p_dbl_Tw+p_dbl_Wr2);
	printf("Lk2=%lf\n",p_dbl_Lk2);
	printf("Writhe1=%lf\n",p_dbl_Wr);
	printf("Tw+Wr1=%lf\n",p_dbl_Tw+p_dbl_Wr);
	printf("Lk=%lf\n",p_dbl_Lk);
#endif
#ifdef PRINT_TO_FILE
	fstream filestr("Mresults7.txt", ios::out|ios::app);
 	filestr<<"Wr_1=\t"<<p_dbl_Wr<<"\tWr_2=\t"<<p_dbl_Wr2<<"\n";
	filestr<<"Wr_1+Tw=\t"<<p_dbl_Wr+p_dbl_Tw<<"\tWr_2+Tw=\t"<<p_dbl_Wr2+p_dbl_Tw<<"\n";
	filestr<<"Lk_1=\t"<<p_dbl_Lk<<"\tLk_2=\t"<<p_dbl_Lk2<<"\n\n\n";
  	filestr.close();
#endif
	
}
void GetFile()
{
	printf("print file name\t");
	scanf("%s",&m_chr_filename);
	FILE *KNOT;
	if((KNOT=fopen( m_chr_filename, "r"))==NULL)
	{
		printf("ERROR coud not open file\n");
		exit(1);
	}
	
	fseek(KNOT, 0L, SEEK_SET);
	int NOL=0;
	int CheckEOF=1;
	fscanf(KNOT,"%d\n",&NOL);		
	fclose(KNOT);
	m_int_N=(int)(NOL);
	m_int_N2=(int)(m_int_N*2);
#ifdef PRINT_TO_FILE
	fstream filestr("Mresults7.txt", ios::out|ios::app);
	filestr<<"Mobius:\tfilename:\t"<<m_chr_filename<<"\n"<<"# of Carbons=\t"<<m_int_N<<"\n";
	filestr.close();
#endif
	//printf("# of Carbons=%d\n",m_int_N);
}
int main()
{	
	GetFile();
	Initialize();
	GetData();
	Topology();
	finish();
	return 1;
}