/***********************************************************/
  /* baryon2000.c                                            */
  /* Copyright (C) 2000 Chiaki Itoh, Meijigakuin University. */
  /* All rights reserved.                                    */
  /***********************************************************/

#include <stdio.h>
#include <math.h>
double sinu[2502],x[21];
double yy[44],ex[44];
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);
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]=5.6369130; /* d-u mass difference */
 x[2]=0.4455472; /* Alpha_s */
 x[3]=4.8772161; /* K*10^-7MeV^3 */
 x[4]=0.2810912; /* Kappa*10^8Mev^-3 */
 x[5]=0.6349005; /* u/s */
 x[6]=1952.9714235;/* c */
 x[7]=435.9089853; /* u */
  x[8]=0.0;
 
 /* experimental values */
 ex[1]=1230.9; /* Delta++ */
 ex[2]=2452.8; /* Sigmac++ */
 ex[3]=2519.4; /* Sigmac*++ */
 ex[4]=938.27231; /* p */
 ex[5]=1231.6; /* Delta+ */
 ex[6]=1189.37;/* Sigma+ */
 ex[7]=1382.8; /* Sigma*+ */
 ex[8]=2453.6; /* Sigmac+ */
 ex[9]=2284.9; /* Lambdac+ */
 ex[10]=939.56563; /* n */
 ex[11]=1233.3; /* Delta0 */
 ex[12]=2452.2; /* Sigmac0 */
 ex[13]=1197.449;/* Sigma- */
 ex[14]=1387.2; /* Sigma*- */
 ex[15]=1192.642;/* Sigma0 */
 ex[16]=1383.7; /* Sigma*0 */
 ex[17]=1115.683;/* Lambda */
 ex[18]=1321.32; /* Xi- */
 ex[19]=1535.0; /* Xi*- */
 ex[20]=1672.45;/* Omega */
 ex[21]=1314.9; /* Xi0 */
 ex[22]=1531.8; /* Xi*0 */
 ex[23]=2573.4; /* Xic+ c{su} */
 ex[24]=2465.6; /* Xic+ c[su] */
 ex[25]=2644.6; /* Xic*+ {csu} */
 ex[26]=2577.3; /* Xic0 c{sd} */
 ex[27]=2470.3; /* Xic0 c[sd] */
 ex[28]=2643.8; /* Xic*0 {csd} */
 ex[29]=2704.0; /* Omegac0 c{ss}*/
 ex[30]=2517.5; /* Sigmac*0 */

 MKAI=0;
  IN=0;
  N=7;
  NK=6;
  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);
  
/*  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(fabs(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((fabs(del[6])>0.1e-13)||(fabs(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\n",IN,x[5],x[6],x[7],FX);
  AV = yy[4]-6.50;
  for(i=1; i<=40; i++){
   yy[i] = yy[i]-AV;
  printf("%13.1f",yy[i]);
  }
  printf("\n");
  for(i=1; i<=30; i++){
   yy[i] = yy[i]+ex[i];
  printf("%13.1f",yy[i]);
  }
  printf("\n");
  for(i=31; i<=40; i++){
   printf("%13.1f",yy[i]);
  }
  printf("\n");
  printf("%13.6f\n",AV);

}

double S(int N)
/*int N;*/
{
  double X,FA,EP,VUD,AMU,AMD,AMS,AMC,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];
/* AMM=1.0+2.0*ALS/(3.0*PI);
 AMMS=1.0-ALS/(12.0*PI);
*/
  AMU=x[7];
  AMS=AMU/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[1]=Z(3)-ex[1];
/* SIGMA C++ */
  am[1]=AMC;
 yy[2]=Z(1)-ex[2];
 yy[3]=Z(3)-ex[3];
/* PROTON  */
  am[1]=AMD;
  q[1]=Q2;
 yy[4]=Z(1)-ex[4];
 yy[5]=Z(3)-ex[5];
/* SIGMA + */
  am[1]=AMS;
  yy[6]=Z(1)-ex[6];
  yy[7]=Z(3)-ex[7];
/* SC + */
  am[1]=AMC;
  am[3]=AMD;
  q[1]=Q1;
  q[3]=Q2;
 yy[8]=Z(1)-ex[8];
  yy[31]=Z(3);
/* LC +  */
  yy[9]=Z(2)-ex[9];
/* NEUTRON  */
  am[1]=AMU;
  am[2]=AMD;
  q[2]=Q2;
 yy[10]=Z(1)-ex[10];
  VUD=W;
 yy[11]=Z(3)-ex[11];
/* SC 0  */
  am[1]=AMC;
 yy[12]=Z(1)-ex[12];
 yy[30]=Z(3)-ex[30];
  am[1]=AMD;
  q[1]=q[2];
 yy[32]=Z(3);
/* SIGMA - */
  am[1]=AMS;
  q[1]=Q2;
 yy[13]=Z(1)-ex[13];
 yy[14]=Z(3)-ex[14];
/* SIGMA 0 */
  am[2]=AMU;
  q[2]=Q1;
 yy[15]=Z(1)-ex[15];
 yy[16]=Z(3)-ex[16];

/* LAMBDA */
 yy[17]=Z(2)-ex[17];
/* XI - */
  am[1]=AMD;
  am[2]=AMS;
  am[3]=AMS;
  q[2]=Q2;
  yy[18]=Z(1)-ex[18];
  yy[19]=Z(3)-ex[19];
/* OMEGA */
  am[1]=AMS;
  yy[20]=Z(3)-ex[20];
/* XI 0  */
  am[1]=AMU;
  q[1]=Q1;
 yy[21]=Z(1)-ex[21];
 yy[22]=Z(3)-ex[22];
/* CSU   */
  am[1]=AMC;
  am[2]=AMS;
  am[3]=AMU;
  q[1]=Q1;
  q[2]=Q2;
  q[3]=Q1;
 yy[23]=Z(1)-ex[23];
 yy[24]=Z(2)-ex[24];
 yy[25]=Z(3)-ex[25];
/*  goto c101; */
/* CSD */
  am[3]=AMD;
  q[3]=Q2;
  yy[26]=Z(1)-ex[26];
  yy[27]=Z(2)-ex[27];
  yy[28]=Z(3)-ex[28];
/* CSS */
  am[3]=AMS;
 yy[29]=Z(1)-ex[29];
 /*goto c110; */
 yy[33]=Z(3);
/* SCC */
 c101:  am[1]=AMS;
  am[2]=AMC;
  am[3]=AMC;
  q[1]=Q2;
  q[2]=Q1;
  q[3]=Q1;
 yy[34]=Z(1);
  yy[35]=Z(3);
/*  goto c102; */
/* DCC   */
  am[1]=AMD;
  yy[36]=Z(1);
  yy[37]=Z(3);
/* UCC   */
  am[1]=AMU;
  q[1]=Q1;
  yy[38]=Z(1);
  yy[39]=Z(3);
/* CCC */
 c102:  am[1]=AMC;
  am[2]=AMC;
  am[3]=AMC;
  q[1]=Q1;
  q[2]=Q1;
  q[3]=Q1;
 yy[40]=Z(3);
/* GOTO 100;   */
c110: FA = 0.0;
  for(i=1; i<=29; i++){
   for(j=i+1; j<=30; j++){
    AV=yy[i]-yy[j];
    FA += AV*AV;
   }
  }

  c100: 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,QQQ,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;
  IMAX=2;
  if(am[2]==am[3])
  goto c35;
  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))));
*/
  goto c38;
 c35:  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];
 c38:  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==3)
  goto c36;
  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];
  goto c37;
 c36:  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];
 c37:  r[1]=rr[1]-RK*dv[1];
  r[2]=rr[2]-RK*dv[2];
 c55:  W=2.0*pow(r[1]/r[2],2.0)-1.0;
  if (IMAX!=2)
  goto c5;
  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];
  goto c7;
   c5: 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;
  c7: 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;
/* printf("%13.6f\n" DMM); */
  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",r[1],r[2],r[3],SZZ,SII);
*/
  return Y;
}

double SUM(double DMM)
/* 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 DMM,XA; */
{
  double XX;
  XX=XA*XA;
  return (XX*exp(-XX)*sqrt(XX+DMM));
}