C语言——一些处理脚本(14)——(待完善)

# include <stdio.h>
# include <math.h>
# include <string.h>
# include <stdlib.h>

# define N   30
//# define NX   500

//double th0=2.182;
//double ph0=1.658;
double l=4.0;

# define L0   14
# define L1   15

int main(void)
{
  int i, j, k, m, n,l;
  double x[N], y[N], z[N];
  double xm1=0, ym1=0,zm1=0;
  double xm2=0, ym2=0,zm2=0;
  double A1,  B1,  C=1;
  double A2,  B2;
  double  AA1=0,BB1=0,CC1=0,DD1=0,AA2=0,BB2=0,CC2=0,DD2=0;
  double L[N], D[N-2], SW1[N], SW2[N], CW1[N], CW2[N], PH1[N], PH2[N], APH1[N], APH2[N],TPH1=0, TPH2=0, TPH3=0;
  char Istring[50];
  FILE *fp;
 
      printf("k=");
      scanf("%d",&k);
      sprintf(Istring,"snp%d",k);
      fp=fopen(Istring,"r");
      for(i=0;i<N;i++)
      {
          fscanf(fp,"%lf %lf %lf",&x[i],&y[i],&z[i]);
      }
      fclose(fp); 

      for(i=0;i<L0;i++)
      {
          xm1+=x[i];
          ym1+=y[i];
          zm1+=z[i];
      }
      xm1=xm1/L0;ym1=ym1/L0;zm1=zm1/L0;

      for(i=L1+1;i<N;i++)
      {
         xm2+=x[i];
         ym2+=y[i];
         zm2+=z[i];
      }
      xm2=xm2/L0;ym2=ym2/L0;zm2=zm2/L0;
  

     for(i=0;i<L0;i++)
    {
        AA1+=(xm1-x[i])*(xm1-x[i]);
        BB1+=(xm1-x[i])*(zm1-z[i]);
        CC1+=(ym1-y[i])*(ym1-y[i]);
        DD1+=(ym1-y[i])*(zm1-z[i]);
    }
    A1=AA1/BB1;B1=CC1/DD1;

    for(i=L1+1;i<N;i++)
    {
        AA2+=(xm2-x[i])*(xm2-x[i]);
        BB2+=(xm2-x[i])*(zm2-z[i]);
        CC2+=(ym2-y[i])*(ym2-y[i]);
        DD2+=(ym2-y[i])*(zm2-z[i]);
    }
    A2=AA2/BB2;B2=CC2/DD2;


   for(i=0;i<N;i++)
   {
      L[i]=sqrt((x[i+1]-x[i])*(x[i+1]-x[i])+(y[i+1]-y[i])*(y[i+1]-y[i])+(z[i+1]-z[i])*(z[i+1]-z[i]));
   }

  for(j=0;j<N-2;j++)
  {
      D[j]=sqrt((x[j+2]-x[j])*(x[j+2]-x[j])+(y[j+2]-y[j])*(y[j+2]-y[j])+(z[j+2]-z[j])*(z[j+2]-z[j]));
  }
  
  for(k=0;k<L0-1;k++)//
  {
      SW1[k]=(fabs(A1*(x[k+1]-x[k])+B1*(y[k+1]-y[k])+C*(z[k+1]-z[k])))/(sqrt(A1*A1+B1*B1+C*C))/L[k];
  }

  for(k=L1+1;k<N-1;k++)//
  {
      SW2[k]=(fabs(A2*(x[k+1]-x[k])+B2*(y[k+1]-y[k])+C*(z[k+1]-z[k])))/(sqrt(A2*A2+B2*B2+C*C))/L[k];
  }

  for(m=0;m<L0-1;m++)//
  {
      CW1[m]=sqrt(1-SW1[m]*SW1[m]);
  }

  for(m=L1+1;m<N-1;m++)//
  {
      CW2[m]=sqrt(1-SW2[m]*SW2[m]);
  }

  for(n=0;n<L0-1;n++)//
  {
      PH1[n]=(L[n]*L[n]+L[n+1]*L[n+1]+2*L[n]*L[n+1]*SW1[n]*SW1[n+1]-D[n]*D[n])/(2*L[n]*L[n+1]*CW1[n]*CW1[n+1]);
      APH1[n]=acos(PH1[n]);
      APH1[n]=acos(-1.)-APH1[n];
  }

  for(n=L1+1;n<N-1;n++)//
  {
      PH2[n]=(L[n]*L[n]+L[n+1]*L[n+1]+2*L[n]*L[n+1]*SW2[n]*SW2[n+1]-D[n]*D[n])/(2*L[n]*L[n+1]*CW2[n]*CW2[n+1]);
      APH2[n]=acos(PH2[n]);
      APH2[n]=acos(-1.)-APH2[n];
  }
  
  for(i=0;i<L0-2;i++)//
  {
  printf("%lf  %lf  %lf  %lf\n",PH1[i],APH1[i],SW1[i],CW1[i]);
  }
  printf("\n");

  for(i=L1+1;i<N-2;i++)//
  {
  printf("%lf  %lf  %lf  %lf\n",PH2[i],APH2[i],SW2[i],CW2[i]);
  }
  printf("\n");

  for(i=0;i<L0-2;i++)
  {
      TPH1+=APH1[i];
  }

  for(i=L1+1;i<N-2;i++)
  {
      TPH2+=APH2[i];
  }

  TPH3=TPH1+TPH2;
  TPH3=TPH3/24;

  printf("TPH3:%lf   %lf   %lf\n",TPH3,TPH3*180/acos(-1.),360/(TPH3*180/acos(-1.))); 
 
  return 0;
}