三瞬子函数集.c

这里面包括了地球物理方面常见的三个程序

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

	int maxx(tb,le,amxx)
	float tb[],*amxx;
	int le;
	{
	  int i;
	  *amxx=tb[1];
	  for(i=1;i<=le;i++)
      if(*amxx<tb[i])  *amxx=tb[i];
	  return;
	}

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

 int    printt(a,l)
 float  a[];
 int    l;
    {
      int i;
      printf("\n output data");
      for(i=1;i<=l;i++)
      {
       printf(" %7.1f", a[i]);
       if(i/5*5==i)  printf("\n");
      }
    }


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

 int    printb(kss,lss,mz)
 int    kss[],lss[];
 int    mz;
   {
    int i;
    for(i=1;i<=mz;i++)
     {
       printf(" %7d", kss[i]);
       printf(" %7d", lss[i]);
       if(i/5*5==i)  printf("\n");
     }
   }


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

 int     phase0(lsm,ksm,ta,le)
 float   lsm[],ksm[],ta[];
 int     le;
     {
       int i;
       float pi,pid,a,b;
       pi=3141.59;
       pid=1570.795;
       for(i=1;i<=le;i++)
	  {
	    a=lsm[i];
	    b=ksm[i];
	    if(a!=0.0)  goto  lp13;
	    if(b<0.0)   goto  lp14;
	    ta[i]=0.0;
	    continue;
lp14:       ta[i]=pi;
	    continue;
lp13:       if(b!=0.0)  goto  lp15;
	    ta[i]=pid;
	    if(a<0.0)  ta[i]=-pid;
	    continue;
lp15:       ta[i]=atan(a/b)*1000.0;
	    if(b<0.0)  ta[i]=ta[i]+pi;
	  }
     }

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

int    freq(ta,tb,le,si)
float  ta[],tb[];
int    le,si;
   {
     int i,jf;
     float dpi,den,fa,fb,dx,dy,dz,x,y,z,del;
     dpi=6283.18;
     den=dpi*si;
     jf=le-1;
     for(i=1;i<=le;i++)
       {
	 fa=ta[i];
	 fb=ta[i+1];
	 dx=fb-fa;
	 dy=fb+dpi-fa;
	 dz=fb-dpi-fa;
	 x=abs(dx);
	 y=abs(dy);
	 z=abs(dz);
	 del=dx;
	 if((fb<fa) && (x>y))  del=dy;
	 if((fb>fa) && (x>z))  del=dz;
	 tb[i]=-del/den*1000000.0;
       }
      tb[le]=tb[jf];
   }

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

 int    hilbt1(hil,nhil)
 float  hil[];
 int    nhil;
    {
	int np,i,nh,icen;
	float  a,a2,cons,an,hh;
	a=36.64677993;
	a2=a*a;
	cons=a2/3.14159265*2.0;
	nh=nhil/2;
	for(np=1;np<=nhil;np++)
	hil[np]=0;
	icen=nh+1;
	an=1.0;
	for(i=1;i<=nh;i+=2)
	 {
	   hh=cons/(an*(an*an+a2));
	   hil[icen+i]=hh;
	   hil[icen-i]=-hh;
	   an=an+2.0;
	 }
    }

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

 int    hilbt2(h,nk)
 float  h[];
 int    nk;
     {
       int np,k1,i2,i1;
       float c,za,zb;
       for(np=1;np<=nk;np++)
       h[np]=0;
       k1=nk;
       i2=k1/2;
       if(k1!=2*i2)  goto  lp415;
       k1++;
       i2=k1/2;
lp415: c=1.0/3.14159236;
       for(i1=1;i1<=i2;i1++)
	 {
	   za=-1;
	   zb=1.0;
	   for(np=1;np<=i1;np++)
	   zb=zb*za;
	   h[i1+i2+1]=c*(1.0-zb)/(float)(i1);
	   h[i2-i1+1]=-h[i2+i1+1];
	  }
     }

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

 int    window(h,nk)
 float  h[];
 int    nk;
    {
       int k,mi;
       float ss;
       for(k=1;k<=nk;k++)
       h[k]=0;
       mi=nk/2+1;
       for(k=1;k<=nk;k++)
	{
	  ss=(2.0*3.14159236*(k-mi))/(float)(nk);
	  h[k]=(0.56+0.44*cos(ss));
	}
    }

/***********************************************/
 int    cvbx(ke,kf,ksm,kc12,nech)
 float  ke[],kf[],ksm[];
 int    kc12,nech;
    {
      int  nh,nchb,np,i,j;
      float sum;
      nh=kc12/2+1;
      nchb=nech-kc12+1;
      for(np=1;np<=nech;np++)
      ksm[np]=0;
      for(i=1;i<=nchb;i++)
       {
	 sum=0.0;
	 for(j=1;j<=kc12;j++)
	 sum+=ke[i+j-1]*kf[j];
	 ksm[i+nh-1]=sum;
       }
    }

/******************************************/
 int    cvbx0(ke,kf,ksm,kc12,nech)
 float  ke[],kf[],ksm[];
 int    kc12,nech;
    {
	int nh,nchb,np,i,j;
	float sum;
	nh=kc12/2+1;
	nchb=nech-kc12+1;
	for(np=1;np<=nech;np++)
	ksm[np]=0;
	for(i=1;i<=nchb;i++)
	 {
	   sum=0;
	   for(j=1;j<=kc12;j++)
	   sum=sum+ke[i+kc12-j]*kf[j];
	   ksm[i+nh-1]=sum;
	 }
    }

/*********************************************/
 int    smooth(ks,ke,nn,nt)
 float  ks[],ke[];
 int    nn,nt;
    {
	int n4,nq,l1,m,k,j;
	float sum;
	n4=nn-nt;
	nq=nt/2+1;
	for(l1=1;l1<=nn;l1++)
	ke[l1]=0;
	for(m=nt;m<=n4;m++)
	 {
	   sum=0;
	   for(k=1;k<=nt;k++)
	   sum=sum+ks[m+k-nq];
	   ke[m]=sum/(float)nt;
	 }
	for(j=1;j<=nn;j++)
	ks[j]=ke[j];
    }

/********************************************/
 int    weight(l,xr,a,ww,mk)
 float  xr[],a[],ww[];
 int    l,mk;
     {
	int ip,kk,j,la,lb,np;
	float e,g;
	if(mk==0)  return;
	for(np=1;np<=l;np++)
	ww[np]=0.0;
	ip=mk/2+1;
	for(j=ip;j<=l;j++)
	 {
	    e=0.0;
	    g=0.0;
	    for(lb=1;lb<=mk;lb++)
	     {
	       la=j+lb-ip;
	       if(la<=0 || la>l)  continue;
	       e=e+a[la]*xr[la];
	       g=g+a[la];
	     }
	    if(g==0.0)  g=1.0;
	    ww[j]=e/g;
	 }
	for(kk=1;kk<=l;kk++)
	xr[kk]=ww[kk];
    }

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

 int    phase1(lsm,ksm,ta,le)
 float  lsm[],ksm[],ta[];
 int    le;
   {
     int i;
     float pii,pidd,a,b,cpk;
     pii=10000.0;
     pidd=9000.0;
     for(i=1;i<=le;i++)
       {
	 a=lsm[i];
	 b=ksm[i];
	 if(a!=0)   goto  lp142;
	 ta[i]=0.0;
	 continue;
lp142:   if(b!=0)   goto  lp144;
	 ta[i]=pidd;
	 if(a<0.0)  ta[i]=-pidd;
	 continue;
lp144:cpk=atan(a/b)*pii/3.141593;
	 cpk=abs(cpk);
	 ta[i]=cpk;
	 if(a<0.0)  ta[i]=-cpk;
       }
     for(i=1;i<=le;i++)
     ta[i]=pidd-ta[i];
   }

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