// CnsoleApplication1.cpp : アプリケーションのエントリ ポイントを定義します。

//



/***********************************************************/

/* baryon2021.cpp                                          */

/* Visual C++:: Console Applecation                        */

/* Copyright (C) 2021 Chiaki Itoh, Meijigakuin University. */

/* All rights reserved.                                    */

/***********************************************************/



//#include <stdio.h>

//#include <math.h>



#include "stdafx.h"

#include <iostream>

#include <cstdio>

using namespace std;



	

	double sinu[2502], x[21];

	double yy[60], ex[60];

	double AV;

	double dv[4], ALS, RC, RK, dom[4], AMM, AMMS, DDD, ALALA, am[4], q[4], W, RO;

	double S(int N);

	double Z(int IZ);

	double SUM(double DMM);

	double F(double DMM, double XA);

	int main()

	{

		double FX, FY, XIN, DL, FZ, DMM;

		double gx[11], del[11];

		int i, N, NK, MKAI, IN, NP, NPP, NPX, NN, II;

		

		DMM = 0.0;

		for (i = 1; i <= 601; i++) {

			DMM += 0.01;

			sinu[i] = SUM(DMM);

		}

		DMM = 10.0;

		for (i = 1001; i <= 2501; i++) {

			DMM += 0.01;

			sinu[i] = SUM(DMM);

		}



		/* input values */

		x[1] = 6.1644051; /* d-u mass difference */

		x[2] = 0.4598088; /* Alpha_s */

		x[3] = 5.2337475; /* K*10^-7MeV^3 */

		x[4] = 0.2104068; /* Kappa*10^8MeV^-3 */

						  /*x[4]=31.4159265358979/(24.0*x[3]);*/

		x[5] = 697.8121636; /* s */

		x[6] = 1969.3908646;/* c */

		x[7] = 443.0130720; /* u */

		x[8] = 5344.0229759; /* b */



							 /* experimental values */

		ex[1] = 938.272029; /* p */

		ex[2] = 939.565360; /* n */

		ex[3] = 1115.683;/* Lambda */

		ex[4] = 1189.37;/* Sigma+ */

		ex[5] = 1192.642;/* Sigma0 */

		ex[6] = 1197.449;/* Sigma- */

		ex[7] = 1314.86; /* Xi0 */

		ex[8] = 1321.71; /* Xi- */

		ex[9] = 1230.55; /* Delta++ */

		ex[10] = 1226.4; /* Delta+ */

		ex[11] = 1233.3; /* Delta0 */

		ex[12] = 1232.6; /* Delta- */

		ex[13] = 1382.8; /* Sigma*+ */

		ex[14] = 1383.7; /* Sigma*0 */

		ex[15] = 1387.2; /* Sigma*- */

		ex[16] = 1531.8; /* Xi*0 */

		ex[17] = 1535.0; /* Xi*- */

		ex[18] = 1672.45;/* Omega */

		ex[19] = 2286.46; /* Lambdac+ */

		ex[20] = 2453.97; /* Sigmac++ */

		ex[21] = 2452.9; /* Sigmac+ */

		ex[22] = 2453.75; /* Sigmac0 */

		ex[23] = 2467.89; /* Xic+ c[su] */

		ex[24] = 2470.99; /* Xic0 c[sd] */

		ex[25] = 2578.2; /* Xic+ c{su} */

		ex[26] = 2578.7; /* Xic0 c{sd} */

		ex[27] = 2695.2; /* Omegac0 c{ss} */

		ex[28] = 2518.41; /* Sigmac*++ */

		ex[29] = 2517.5; /* Sigmac*+ {cdu} */

		ex[30] = 2518.48; /* Sigmac*0 */

		ex[31] = 2645.10; /* Xic*+ {csu} */

		ex[32] = 2646.16; /* Xic*0 {csd} */

		ex[33] = 2765.9; /* Omegac0 {css} */

		ex[34] = 5619.60;/* Lambdab0 */

		ex[35] = 5810.56; /* Sigmab+ */

		ex[36] = 5809.9; /* Sigmab0 */

		ex[37] = 5815.64; /* Sigmab- */

		ex[38] = 5791.9; /* Xib0 */

		ex[39] = 5797.0; /* Xib- */

		ex[40] = 5924.2; /* Xib0 */

		ex[41] = 5935.02; /* Xib- */

		ex[42] = 6030.7; /* Omegab- */

		ex[43] = 5830.32; /* Sigmab*+ */

		ex[44] = 5839.9; /* Sigmab*0 */

		ex[45] = 5835.02; /* Sigmab*- */

		ex[46] = 5952.3; /* Xib*0 */

		ex[47] = 5955.33; /* Xib*- */

		ex[48] = 6059.3; /* Omegab- */



		MKAI = 0;

		IN = 0;

		N = 8;

		NK = 7;

		FX = S(N);



		printf("%d %13.6f %13.6f %13.6f %13.6f\n",IN,x[1],x[2],x[3],x[4]);

		printf("%d %13.6f %13.6f %13.6f %13.6f\n",IN,x[5],x[6],x[7],FX);

		/*

		 FZ=FX+100;

		goto c81;

		*/

		FZ = FX;

	c44: DL = 0.01;

	c45: IN = 0;

		NP = 1;

		NPP = 0;

		NPX = N + 1;

		for (i = 1; i <= N; i++)

			del[i] = 1.e-4;

	c43: NN = 0;

		for (i = 1; i <= N; i++)

			del[i] *= 10.0;

		if (IN != 0)

			goto c72;

		FX = S(N);

		goto c72;

	c50: FX = FY;

		del[IN] *= 1.2;

	c72:  IN++;

		if (IN == NPX)

			IN = 1;

		if (IN == NP)

			goto c61;

	c73: XIN = x[IN];

	c75: x[IN] = (1.0 + del[IN])*XIN;

		FY = S(N);

		if (FX >= FY)

			goto c50;

		del[IN] *= -0.5;

		if (abs(del[IN])>0.1e-13)

			goto c75;

		x[IN] = XIN;

		goto c72;

	c61: if (NPP != 0)

		goto c62;

		 NPP = 1;

		 goto c73;

	 c62: for (i = 1; i <= N; i++) {

		 gx[i] = x[i];

		 x[i] *= (1.0 + del[i]);

	 }

		  FY = S(N);

		  if (FX<FY)

			  goto c41;

	  c15: for (i = 1; i <= N; i++)

		  del[i] *= 1.1;

		   FX = FY;

		   goto c62;

	   c41: for (i = 1; i <= N; i++)

		   x[i] = gx[i];

			del[IN] = -del[IN];

			x[IN] *= (1.0 + del[IN]);

			FY = S(N);

			if (FX <= FY)

				goto c16;

			/* printf("%d %13.6f %13.6f %13.6f %13.6f\n",IN,x[1],x[2],x[3],x[4]);

			printf("%d %13.6f %13.6f %13.6f %13.6f\n",IN,x[5],x[6],x[7],FY);

			*/

			goto c15;

		c16: x[IN] = gx[IN];

			del[IN] = -del[IN];

			if (NN<10)

				goto c51;

			/* printf("%d %13.6f %13.6f %13.6f %13.6f\n",IN,x[1],x[2],x[3],x[4]);

			printf("%d %13.6f %13.6f %13.6f %13.6f\n",IN,x[5],x[6],x[7],FX);

			*/

			goto c43;

		c51: NN++;

			NPP = 0;

			NP++;

			if (NP == NPX)

				NP = 1;

			if ((abs(del[6])>0.1e-13) || (abs(del[4])>0.1e-13))

				goto c72;

			if (FX<FZ)

				goto c80;

			for (i = 1; i <= N; i++) {

				II = i + N;

				x[i] = x[II];

			}

			if (DL>0.0)

				goto c83;

			DL *= -0.5;

			if (DL>1.e-3)

				goto c85;

			NK--;

			if (NK == 0)

			{

				MKAI++;

				if (MKAI = 10)

					goto c81;

				NK = 6;

				goto c44;

			}

			goto c44;

		c83: DL = -DL;

			goto c85;

		c80:  printf("%d %13.6f %13.6f %13.6f %13.6f\n",NK,x[1],x[2],x[3],x[4]);

			 printf("%d %13.6f %13.6f %13.6f %13.6f\n",IN,x[5],x[6],x[7],FX);



			FZ = FX;

			for (i = 1; i <= N; i++) {

				II = i + N;

				x[II] = x[i];

			}

		c85: x[NK] *= (1.0 + DL);

			goto c45;

		c81:

			if (FX>FZ)

			{

				FX = FZ;

				for (i = 1; i <= N; i++) {

					II = i + N;

					x[i] = x[II];

				}

				FX = S(N);

			}

	        printf("%d %13.7f %13.7f %13.7f %13.7f\n",NK,x[1],x[2],x[3],x[4]);

				printf("%d %13.7f %13.7f %13.7f %13.7f %13.7f\n",IN,x[5],x[6],x[7],x[8],FX);



			AV = yy[1] - 5.24;



			for (i = 1; i <= 8; i++) {

				yy[i] = yy[i] - AV;

				printf("%13.1f",yy[i]);

			}

			printf("\n");



			for (i = 1; i <= 8; i++) {

				yy[i] = yy[i] + ex[i];

				printf("%13.1f",yy[i]);

			}

			printf("\n");



			for (i = 9; i <= 16; i++) {

				yy[i] = yy[i] - AV;

				printf("%13.1f",yy[i]);

			}

			printf("\n");



			for (i = 9; i <= 16; i++) {

				yy[i] = yy[i] + ex[i];

				printf("%13.1f",yy[i]);

			}

			printf("\n");



			for (i = 17; i <= 24; i++) {

				yy[i] = yy[i] - AV;

				printf("%13.1f",yy[i]);

			}

			printf("\n");



			for (i = 17; i <= 24; i++) {

				yy[i] = yy[i] + ex[i];

				printf("%13.1f",yy[i]);

			}

			printf("\n");

			

			for (i = 25; i <= 32; i++) {

				yy[i] = yy[i] - AV;

				printf("%13.1f",yy[i]);

			}

			printf("\n");



			for (i = 25; i <= 32; i++) {

				yy[i] = yy[i] + ex[i];

				printf("%13.1f",yy[i]);

			}

			printf("\n");



			for (i = 33; i <= 40; i++) {

				yy[i] = yy[i] - AV;

				printf("%13.1f",yy[i]);

			}

			printf("\n");



			for (i = 33; i <= 40; i++) {

				yy[i] = yy[i] + ex[i];

				printf("%13.1f",yy[i]);

			}

			printf("\n");



			for (i = 41; i <= 48; i++) {

				yy[i] = yy[i] - AV;

				printf("%13.1f",yy[i]);

			}

			printf("\n");



			for (i = 41; i <= 48; i++) {

				yy[i] = yy[i] + ex[i];

				printf("%13.1f",yy[i]);

			}

			printf("\n");



			printf("\n");

			printf("%13.2f %13.2f\n",AV,yy[17]-yy[8]-(yy[15]-yy[6]));

			

		/*	for (i = 50; i <= 56; i++) {

				yy[i] = yy[i] -AV;

				printf("%13.1f", yy[i]);

			}

			printf("\n");

        */

			return 0;

	}



	double S(int N)

		/*int N;*/

	{

		double FA, EP, VUD, AMU, AMD, AMS, AMC, AMB, Q1, Q2, PI;

		int i, j;

		PI = 3.141592653589793;

		EP = x[1];

		ALS = x[2];

		RO = x[3] * 1.e7;

		RK = x[4] / 1.e8;

		AMC = x[6];

		AMB = x[8];

		/* AMM=1.0+2.0*ALS/(3.0*PI);

		AMMS=1.0-ALS/(12.0*PI);

		*/

		AMU = x[7];

		AMS = x[5];

		/*RS=x[8]*1.e6;*/

		AMD = AMU + EP;

		/* ALALA=2.0*ALS*ALS*(1.0/(AMU*AMU)+1.0/(AMD*AMD)+1.0/(AMS*AMS));

		*/

		Q1 = 2.0 / 3.0;

		Q2 = -1.0 / 3.0;

		/* DELTA ++ */

		am[1] = AMU;

		am[2] = AMU;

		am[3] = AMU;

		q[1] = Q1;

		q[2] = Q1;

		q[3] = Q1;

		yy[9] = Z(3) - ex[9];

		/* SIGMA C++ */

		am[1] = AMC;

		yy[20] = Z(1) - ex[20];

		yy[28] = Z(3) - ex[28];

		/* PROTON  */

		am[1] = AMD;

		q[1] = Q2;

		yy[1] = Z(1) - ex[1];

		yy[10] = Z(3) - ex[10];

		/* SIGMA + */

		am[1] = AMS;

		yy[4] = Z(1) - ex[4];

		yy[13] = Z(3) - ex[13];

		/* SC + */

		am[1] = AMC;

		am[3] = AMD;

		q[1] = Q1;

		q[3] = Q2;

		yy[21] = Z(1) - ex[21];

		yy[29] = Z(3) - ex[29];

		/* LC +  */

		yy[19] = Z(2) - ex[19];

		/* NEUTRON  */

		am[1] = AMU;

		am[2] = AMD;

		q[2] = Q2;

		yy[2] = Z(1) - ex[2];

		VUD = W;

		yy[11] = Z(3) - ex[11];

		/* SC 0  */

		am[1] = AMC;

		yy[22] = Z(1) - ex[22];

		yy[30] = Z(3) - ex[30];

		am[1] = AMD;

		q[1] = q[2];

		yy[12] = Z(3) - ex[12];

		/* SIGMA - */

		am[1] = AMS;

		q[1] = Q2;

		yy[6] = Z(1) - ex[6];

		yy[15] = Z(3) - ex[15];

		/* SIGMA 0 */

		am[2] = AMU;

		q[2] = Q1;

		yy[5] = Z(1) - ex[5];

		yy[14] = Z(3) - ex[14];



		/* LAMBDA */

		yy[3] = Z(2) - ex[3];

		/* XI - */

		am[1] = AMD;

		am[2] = AMS;

		am[3] = AMS;

		q[2] = Q2;

		yy[8] = Z(1) - ex[8];

		yy[17] = Z(3) - ex[17];

		/* OMEGA */

		am[1] = AMS;

		yy[18] = Z(3) - ex[18];

		/* XI 0  */

		am[1] = AMU;

		q[1] = Q1;

		yy[7] = Z(1) - ex[7];

		yy[16] = Z(3) - ex[16];

		/* CSU   */

		am[1] = AMC;

		am[2] = AMS;

		am[3] = AMU;

		q[1] = Q1;

		q[2] = Q2;

		q[3] = Q1;

		yy[25] = Z(1) - ex[25];

		yy[23] = Z(2) - ex[23];

		yy[31] = Z(3) - ex[31];

		/* CSD */

		am[3] = AMD;

		q[3] = Q2;

		yy[26] = Z(1) - ex[26];

		yy[24] = Z(2) - ex[24];

		yy[32] = Z(3) - ex[32];

		/* CSS */

		am[3] = AMS;

		yy[27] = Z(1) - ex[27];

		yy[33] = Z(3) - ex[33];

		/* LAMBDAB 0 */

		am[1] = AMB;

		am[2] = AMU;

		am[3] = AMD;

		q[1] = Q2;

		q[2] = Q1;

		q[3] = Q2;

		yy[34] = Z(2) - ex[34];

		yy[36] = Z(1) - ex[36];

		yy[44] = Z(3) - ex[44];

		/* SIGMAB - */

		am[2] = AMD;

		q[2] = Q2;

		yy[37] = Z(1) - ex[37];

		yy[45] = Z(3) - ex[45];

		/* SIGMAB + */

		am[3] = AMU;

		q[3] = Q1;

		am[2] = AMU;

		q[2] = Q1;

		yy[35] = Z(1) - ex[35];

		yy[43] = Z(3) - ex[43];

		/* Xib0 */

		am[2] = AMS;

		q[2] = Q2;

		yy[40] = Z(1) - ex[40];

		yy[38] = Z(2) - ex[38];

		yy[46] = Z(3) - ex[46];

		/* Xib- */

		am[3] = AMD;

		q[3] = Q2;

		yy[41] = Z(1) - ex[41];

		yy[39] = Z(2) - ex[39];

		yy[47] = Z(3) - ex[47];

		/* Omegab- */

		am[3] = AMS;

		yy[42] = Z(1) - ex[42];

		yy[48] = Z(3) - ex[48];



		goto c110;

		RK *= 2.0;

		/* SCC */

		am[1] = AMS;

		am[2] = AMC;

		am[3] = AMC;

		q[1] = Q2;

		q[2] = Q1;

		q[3] = Q1;

		yy[52] = Z(1);

		yy[55] = Z(3);

		/* DCC   */

		am[1] = AMD;

		yy[51] = Z(1);

		yy[54] = Z(3);

		/* UCC   */

		am[1] = AMU;

		q[1] = Q1;

		yy[50] = Z(1);

		yy[53] = Z(3);

		/* CCC */

		am[1] = AMC;

		am[2] = AMC;

		am[3] = AMC;

		q[1] = Q1;

		q[2] = Q1;

		q[3] = Q1;

		yy[56] = Z(3);

	c110: FA = 0.0;

		for (i = 1; i <= 46; i++) {

			for (j = i + 1; j <= 47; j++) {

				if (i != 10 && j != 10 && i != 12 && j != 12 && i != 36 && j != 36 && i != 40 && j != 40 && i != 42 && j != 42 && i != 44 && j != 44) {

					AV = yy[i] - yy[j];

					FA += AV * AV;

					if (i == 19 || j == 19)

						FA += 2.0 * AV * AV;

					if (i == 18 || j == 18 || i == 34 || j == 34)

						FA += AV * AV;

					if (i == 8 || j == 8 || i == 7 || j == 7)

						FA += 0.5*AV*AV;

				}

			}

		}



		return FA;

	}



	double Z(int IZ)

		/* int IZ; */

	{

		double AM2, RINV, SZ, SISJ, RC2, RC4, D, SI, DMM, SR, SD, SSD;

		double A1, A2, A3, A5, A6, PI, C, CK, CM, CM2, CMM, CMS, CQ3;

		double CHM, vv[4], ch[4], rr[4], SUN, SZZ, SII, SDMM, Y;

		double r[4], a[4];

		int i, IMAX, J, K, M;

		PI = 3.141592653589793;

		Y = 0.0;

		A5 = PI / 4.0;



		if (am[2] == am[3]) {

			IMAX = 2;

			AM2 = am[2] + am[3];

			rr[2] = pow(PI*PI / (12.0*RO*am[2]), 0.25);

			A1 = sqrt((am[1] + AM2) / (2.0*am[1] * RO*AM2));

			rr[1] = 0.5*sqrt(PI*(A1 + sqrt(1.0 / (12.0*RO*am[2]))));

			rr[3] = rr[1];

		}

		else {

			IMAX = 3;

			CM = am[1] + am[2] + am[3];

			CM2 = am[1] * am[1] + am[2] * am[2] + am[3] * am[3];

			CMS = am[1] * am[2] * am[3];

			CMM = CM * CM - CM2;

			CK = sqrt(0.5*RO*CMS) / (4.0*CM);

			CQ3 = sqrt(2.0*(CM2*CM2 - pow(am[1], 4.0) - pow(am[2], 4.0) - pow(am[3], 4.0) - 2.0*CMS*CM));

			C = 2.0*(am[1] * am[2] - am[3] * am[3]) / CQ3;

			A1 = CK * (sqrt(CMM - CQ3)*(1.0 + C) + sqrt(CMM + CQ3)*(1.0 - C)) / am[3];

			C = 2.0*(am[2] * am[3] - am[1] * am[1]) / CQ3;

			A2 = CK * (sqrt(CMM - CQ3)*(1.0 + C) + sqrt(CMM + CQ3)*(1.0 - C)) / am[1];

			C = 2.0*(am[3] * am[1] - am[2] * am[2]) / CQ3;

			A3 = CK * (sqrt(CMM - CQ3)*(1.0 + C) + sqrt(CMM + CQ3)*(1.0 - C)) / am[2];

			C = PI / (8.0*(A1*A2 + A2 * A3 + A3 * A1));

			rr[1] = sqrt((A2 + A3)*C);

			rr[2] = sqrt((A3 + A1)*C);

			rr[3] = sqrt((A1 + A2)*C);

			/* AM2=am[2]+am[3];

			A3=am[2]*am[2]+am[2]*am[3]+am[3]*am[3];

			rr[2]=pow(PI*PI*pow(AM2,3.0)/(32.0*RO*am[2]*am[3]*A3),0.25);

			A1=sqrt((am[1]+AM2)/(2.0*am[1]*RO*AM2));

			A2=2.0*RO*AM2*A3;

			rr[1]=0.5*sqrt(PI*(A1+sqrt(pow(am[3],3.0)/(am[2]*A2))));

			rr[3]=0.5*sqrt(PI*(A1+sqrt(pow(am[2],3.0)/(am[3]*A2))));

			*/

		}



		for (i = 1; i <= IMAX; i++) {

			J = (i % 3) + 1;

			K = ((i + 1) % 3) + 1;

			switch (IZ) {

			case 1:

				if (i == 2)

					SISJ = 0.25;

				else

					SISJ = -0.5;

				break;

			case 2:

				if (i == 2)

					SISJ = -0.75;

				else

					SISJ = 0.0;

				break;

			case 3:

				SISJ = 0.25;

			}

			RC = rr[i];

			RC2 = RC * RC;

			RC4 = pow(RC, 4.0);

			ch[i] = -2.0*ALS / 3.0 + q[i] * q[J] / 137.036;

			/* CHM= -2.0*ALS*AMMS*AMMS/3.0+q[i]*q[J]*AMM*AMM/137.036;

			*/

			vv[i] = (1.0 / (am[i] * am[i]) + 1.0 / (am[J] * am[J]))*ch[i] + 16.0*SISJ*ch[i] / (3.0*am[i] * am[J]);

			dv[i] = -ch[i] / RC2 + 3.0*PI*vv[i] / (16.0*RC4);

			/* +ALALA/(6.D1*RC4); */

		}

		if (IMAX == 2) {

			ch[3] = ch[1];

			dv[3] = dv[1];

			SD = rr[1] / rr[2];

			SSD = SD * SD;

			SR = PI / pow(4.0*rr[1] * rr[1] - rr[2] * rr[2], 2.0);

			dv[1] = dv[1] + (-ch[1] * (4.0*SSD + 3.0 - 0.5 / SSD) / (am[1] * am[2]) + ch[2] * SD / (am[2] * am[3]) + ch[3] * (4.0*SSD - 3.0) / (am[3] * am[1]))*SR;

			dv[2] = dv[2] + (ch[2] * (-12.0*SSD*SSD + 7.0*SSD - 1.5) / (am[2] * am[3]) + 2.0*ch[1] / (am[1] * am[2] * SD))*SR;

			dv[3] = dv[1];

		}

		else {

			A6 = 2.0*(pow(rr[1] * rr[2], 2.0) + pow(rr[2] * rr[3], 2.0) + pow(rr[3] * rr[1], 2.0));

			A6 = A6 - pow(rr[1], 4.0) - pow(rr[2], 4.0) - pow(rr[3], 4.0);

			A1 = -rr[1] * rr[1] + rr[2] * rr[2] + rr[3] * rr[3];

			A2 = rr[1] * rr[1] - rr[2] * rr[2] + rr[3] * rr[3];

			A3 = rr[1] * rr[1] + rr[2] * rr[2] - rr[3] * rr[3];

			dv[1] -= ch[1] * PI*A1 / (2.0*am[1] * am[2] * rr[1] * rr[1] * A6);

			dv[1] += ch[1] * PI*(-1.0 - 2.0*A1*A1 / A6) / (am[1] * am[2] * A6);

			dv[1] += ch[2] * PI*((1.0 - 2.0*A1*A2 / A6) / (am[2] * am[3] * A6))*rr[1] / rr[2];

			dv[1] += ch[3] * PI*((1.0 - 2.0*A1*A3 / A6) / (am[3] * am[1] * A6))*rr[1] / rr[3];

			dv[2] -= ch[2] * PI*A2 / (2.0*am[2] * am[3] * rr[2] * rr[2] * A6);

			dv[2] += ch[2] * PI*(-1.0 - 2.0*A2*A2 / A6) / (am[2] * am[3] * A6);

			dv[2] += ch[3] * PI*((1.0 - 2.0*A2*A3 / A6) / (am[3] * am[1] * A6))*rr[2] / rr[3];

			dv[2] += ch[1] * PI*((1.0 - 2.0*A2*A1 / A6) / (am[1] * am[2] * A6))*rr[2] / rr[1];

			dv[3] -= ch[3] * PI*A3 / (2.0*am[3] * am[1] * rr[3] * rr[3] * A6);

			dv[3] += ch[3] * PI*(-1.0 - 2.0*A3*A3 / A6) / (am[3] * am[1] * A6);

			dv[3] += ch[1] * PI*((1.0 - 2.0*A3*A1 / A6) / (am[1] * am[2] * A6))*rr[3] / rr[1];

			dv[3] += ch[2] * PI*((1.0 - 2.0*A3*A2 / A6) / (am[2] * am[3] * A6))*rr[3] / rr[2];

			r[3] = rr[3] - RK * dv[3];

		}

		r[1] = rr[1] - RK * dv[1];

		r[2] = rr[2] - RK * dv[2];

		W = 2.0*pow(r[1] / r[2], 2.0) - 1.0;

		if (IMAX == 2) {

			r[3] = r[1];

			vv[3] = vv[1];

			a[1] = A5 / ((1.0 + 2.0*W)*r[2] * r[2]);

			a[2] = W * a[1];

			a[3] = a[1];

		}

		else {

			A1 = r[1] * r[1] / A5;

			A2 = r[2] * r[2] / A5;

			A3 = r[3] * r[3] / A5;

			A6 = 2.0*(A1*A2 + A2 * A3 + A3 * A1) - A1 * A1 - A2 * A2 - A3 * A3;

			a[1] = (-A1 + A2 + A3) / A6;

			a[2] = (A1 - A2 + A3) / A6;

			a[3] = (A1 + A2 - A3) / A6;

		}

		  SZZ = 0.0;

		  SII = 0.0;

		  for (i = 1; i <= IMAX; i++) {

			  J = (i % 3) + 1;

			  K = ((i + 1) % 3) + 1;

			  RINV = 1.0 / r[i];

			  SZ = (RINV + RINV * 2.0*a[i] / (am[i] * am[J]))*ch[i] - PI * pow(RINV, 3.0)*vv[i] / 16.0 + 0.954929658551372*RO*r[i] * r[i];

			  /* /sqrt(1.0+RS*r[i]*r[i]);

			  +0.5*pow(PI*r[i]/(4.0*RO*RO),2.0);

			  -ALALA/(1.8e2*pow(r[i],3.0));

			  */

			  if ((IMAX == 2) && (i == 1))

				  SZ = 2.0*SZ;

			  D = 0.5 / (a[i] + a[K]);

			  Y += SZ;

			  SZZ += SZ;

			  DMM = am[i] * am[i] * D;

			  SDMM = DMM * 100.0;

			  if ((DMM >= 0.01) && (DMM <= 6.0))

				  goto c63;

			  if ((DMM >= 10.01) && (DMM <= 25.0))

				  goto c63;

				  

			  SUN = SUM(DMM);

			  goto c66;

		  c63: M = floor(SDMM);

			  SUN = (DMM - M / 100.)*(sinu[M + 1] - sinu[M]) / 0.01 + sinu[M];



		  c66: SI = 2.256758334*SUN / sqrt(D);

			  if ((IMAX == 2) && (i == 2))

				  SI *= 2.0;

			  Y += SI;

			  SII += SI;

			  /* printf("%13.6f %13.6f %13.6f %13.6f %13.6f\n",SZ,SI,SUN,DMM,sinu[M]);

			  printf("%13.6f %13.6f\n",sinu[M+1],float(M);

			  */

		  }

		  /* printf("%13.6f %13.6f %13.6f %13.6f %13.6f\n",am[1],am[2],am[3],q[1],q[2]);

		  printf("%13.6f\n",q[3]);

		  printf("%13.6f %13.6f %13.6f %13.6f %13.6f\n",a[1],a[2],a[3],SZZ,SII);

		  */

		  return Y;

	}



	double SUM(double DMM)

	{

		double AI, BI, H, XA, SA;

		int i, NL;

		AI = 0.0;

		BI = 6.0;

		NL = 16;

		H = (BI - AI) / 32.0;

		XA = AI;

		SA = F(DMM, XA) / 2.0;

		for (i = 1; i <= NL; i++) {

			XA += H;

			SA += F(DMM, XA)*2.0;

			XA += H;

			SA += F(DMM, XA);

		}

		SA *= H * 2.0 / 3.0;

		return SA;

	}



	double F(double DMM, double XA)

	{

		double XX;

		XX = XA * XA;

		return (XX*exp(-XX)*sqrt(XX + DMM));

	}