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

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

# define N   30

//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;
  double x[N], y[N], z[N];
  double xm1, ym1,zm1;
  double xm2, ym2,zm2;
  double A1,  B1,  C=1;
  double A2,  B2;
  double  AA1,BB1,CC1,DD1,AA2,BB2,CC2,DD2;
  double L[7], D[6], SW[7], CW[7], PH[6],APH[6], TPH1=0, TPH2=0, TPH3=0;
  FILE *fp;
 
  fp=fopen("snp2","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<N-1;k++)
  {
      SW[k]=(fabs(A*(x[k+1]-x[k])+B*(y[k+1]-y[k])+C*(z[k+1]-z[k])))/(sqrt(A*A+B*B+C*C))/L[k];
  }

  for(m=0;m<N-1;m++)
  {
      CW[m]=sqrt(1-SW[m]*SW[m]);
  }

  for(n=0;n<N-2;n++)
  {
      PH[n]=(L[n]*L[n]+L[n+1]*L[n+1]+2*L[n]*L[n+1]*SW[n]*SW[n+1]-D[n]*D[n])/(2*L[n]*L[n+1]*CW[n]*CW[n+1]);
      APH[n]=acos(PH[n]);
      APH[n]=acos(-1.)-APH[n];
  }
  
  for(i=0;i<L0-2;i++)
  {
      TPH1+=APH[i];
  }

  for(i=L1+1;i<N-2;i++)
  {
      TPH2+=APH[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;
}