打印时遇到fail，求助。

[vfeeok1]

就是出现如图上的问题，打印时出现fail。
但是我调试时，发现循环中每次计算结果均正确。





**************************************************************************************************************
#include "stdio.h"
#include "math.h"
#include "stdlib.h"
#define N 4 //the NO. of the layers
//******************************
float pmax(float *P,int m)
{float max=*P;
int i;
for(i=0;i<m-1;i++)
{if(*(P+i)>max)
max=*(P+i);
}
return max;
}

double fp(float *h,float *P,float x,int m,double p)
{int i;
double f;
f=-0.5*x/p;
for(i=0;i<m-1;i++)
{f=f+(*(h+i))*(*(P+i))/sqrt(1-p*p*(*(P+i))*(*(P+i)));
}
return f;
}

double P_p(float *h, float *P,float x,int m,double eps)
{double p=0.5/(*P),p1=0,p2=1/(*P),f;
float maxp;
maxp=pmax(P,m);
p=0.5/maxp;p2=1/maxp;
f=fp(h,P,x,m,p);
while(fabs(f)>eps)
{if(f>eps)
{p2=p;
p=(p+p1)/2.0;
f=fp(h,P,x,m,p);
}
else
{p1=p;
p=(p+p2)/2.0;
f=fp(h,P,x,m,p);
}
}
return (p);
}

float psmax(float *P,float *S,int m)
{float max=*P;
int i;
for(i=0;i<m-1;i++)
{if(*(P+i)>max)
max=*(P+i);
if(*(S+i)>max)
max=*(S+i);
}
return max;
}

double fs(float *h,float *P,float *S,float x,int m,double p)
{int i;
double f;
f=-x/p;
for(i=0;i<m-1;i++)
{f=f+(*(h+i))*(*(P+i))/sqrt(1-p*p*(*(P+i))*(*(P+i)))+(*(h+i))*(*(S+i))/sqrt(1-p*p*(*(S+i))*(*(S+i)));
}
return f;
}

double S_p(float *h,float *P,float *S,float x,int m,double eps)
{double p,p1=0,p2,f;
float maxps;
maxps=psmax(P,S,m);
p=0.5/maxps;
p2=1/maxps;
f=fs(h,P,S,x,m,p);
while(fabs(f)>eps)
{if(f>eps)
{p2=p;
p=(p+p1)/2.0;
f=fs(h,P,S,x,m,p);
}
else
{p1=p;
p=(p+p2)/2.0;
f=fs(h,P,S,x,m,p);
}
}
return (p);
}

//***********************************
int gaus(a,b,n)
int n;
double a[],b[];
{ int *js,l,k,i,j,is,p,q;
double d,t;
js=malloc(n*sizeof(int));
l=1;
for (k=0;k<=n-2;k++)
{ d=0.0;
for (i=k;i<=n-1;i++)
for (j=k;j<=n-1;j++)
{t=fabs(a[i*n+j]);
if (t>d) { d=t; js[k]=j; is=i;}
}
if (d+1.0==1.0) l=0;
else
{ if (js[k]!=k)
for (i=0;i<=n-1;i++)
{ p=i*n+k; q=i*n+js[k];
t=a[p]; a[p]=a[q]; a[q]=t;
}
if (is!=k)
{ for (j=k;j<=n-1;j++)
{ p=k*n+j; q=is*n+j;
t=a[p]; a[p]=a[q]; a[q]=t;
}
t=b[k]; b[k]=b[is]; b[is]=t;
}
}
if (l==0)
{ free(js); printf("fail\n");
return(0);
}
d=a[k*n+k];
for (j=k+1;j<=n-1;j++)
{ p=k*n+j; a[p]=a[p]/d;}
b[k]=b[k]/d;
for (i=k+1;i<=n-1;i++)
{ for (j=k+1;j<=n-1;j++)
{ p=i*n+j;
a[p]=a[p]-a[i*n+k]*a[k*n+j];
}
b[i]=b[i]-a[i*n+k]*b[k];
}
}
d=a[(n-1)*n+n-1];
if (fabs(d)+1.0==1.0)
{ free(js); printf("fail\n");
return(0);
}
b[n-1]=b[n-1]/d;
for (i=n-2;i>=0;i--)
{ t=0.0;
for (j=i+1;j<=n-1;j++)
t=t+a[i*n+j]*b[j];
b[i]=b[i]-t;
}
js[n-1]=n-1;
for (k=n-1;k>=0;k--)
if (js[k]!=k)
{ t=b[k]; b[k]=b[js[k]]; b[js[k]]=t;}
free(js);
return(1);
}

double *zoep_P(float Pin,float Pout,float Sin,float Sout,float Roin,float Roout,double p,double *temp)
{double A[4][4],B[4];
int i;
A[0][0]=Pin*p;A[0][1]=sqrt(1-p*p*Sin*Sin);A[0][2]=-Pout*p;A[0][3]=sqrt(1-p*p*Sout*Sout);
A[1][0]=-sqrt(1-p*p*Pin*Pin);A[1][1]=Sin*p;A[1][2]=-sqrt(1-p*p*Pout*Pout);A[1][3]=-Sout*p;
A[2][0]=2*Sin*Sin*p*sqrt(1-p*p*Pin*Pin);A[2][1]=Sin*(1-2*Sin*Sin*p*p);A[2][2]=2*Roout/Roin*Sout*Sout*p*sqrt(1-p*p*Pout*Pout);A[2][3]=-Roout/Roin*Sout*(1-2*p*p*Sout*Sout);
A[3][0]=Pin*(1-2*p*p*Sin*Sin);A[3][1]=-2*Sin*Sin*p*sqrt(1-p*p*Sin*Sin);A[3][2]=-Roout/Roin*Pout*(1-2*p*p*Sout*Sout);A[3][3]=-2*Roout/Roin*Sout*Sout*p*sqrt(1-p*p*Sout*Sout);
B[0]=-Pin*p;B[1]=-sqrt(1-p*p*Pin*Pin);B[2]=2*Sin*Sin*p*sqrt(1-p*p*Pin*Pin);B[3]=-Pin*(1-2*p*p*Sin*Sin);
gaus(A,B,4);
for(i=0;i<4;i++)
{*(temp+i)=B[i];}
return temp;
}


double *zoep_S(float Pin,float Pout,float Sin,float Sout,float Roin,float Roout,double p,double *temp)
{double A[4][4],B[4];
int i;
A[0][0]=Pin*p;A[0][1]=sqrt(1-p*p*Sin*Sin);A[0][2]=-Pout*p;A[0][3]=sqrt(1-p*p*Sout*Sout);
A[1][0]=-sqrt(1-p*p*Pin*Pin);A[1][1]=Sin*p;A[1][2]=-sqrt(1-p*p*Pout*Pout);A[1][3]=-Sout*p;
A[2][0]=2*Sin*Sin*p*sqrt(1-p*p*Pin*Pin);A[2][1]=Sin*(1-2*Sin*Sin*p*p);A[2][2]=2*Roout/Roin*Sout*Sout*p*sqrt(1-p*p*Pout*Pout);A[2][3]=-Roout/Roin*Sout*(1-2*p*p*Sout*Sout);
A[3][0]=Pin*(1-2*p*p*Sin*Sin);A[3][1]=-2*Sin*Sin*p*sqrt(1-p*p*Sin*Sin);A[3][2]=-Roout/Roin*Pout*(1-2*p*p*Sout*Sout);A[3][3]=-2*Roout/Roin*Sout*Sout*p*sqrt(1-p*p*Sout*Sout);
B[0]=sqrt(1-p*p*Sin*Sin);B[1]=-Sin*p;B[2]=-Sin*(1-2*Sin*Sin*p*p);B[3]=-2*Sin*Sin*p*sqrt(1-p*p*Sin*Sin);
gaus(A,B,4);
for(i=0;i<4;i++)
{*(temp+i)=B[i];}
return temp;
}

//*****************************************
void P_1x_tpara(float *P,float *S,float *Ro,float *h,float x,int m,double *P_time,double *P_para,double *P_thita) //*
{double p,tP=1.0,*temp,Ptime=0;
int i,it;
double pthita; //*
temp=calloc(4,sizeof(double));
p=P_p(h,P,x,m,0.000000001);
pthita=asin(p*(*P)); //*
for(it=0;it<m-1;it++)
{Ptime=Ptime+2*(*(h+it))/sqrt(1-p*p*(*(P+it))*(*(P+it)))/(*(P+it));
}
if(m==2)
{zoep_P(*P,*(P+1),*S,*(S+1),*Ro,*(Ro+1),p,temp);
tP=*temp;
}
else
{ for(i=0;i<m-2;i++)
{zoep_P(*(P+i),*(P+i+1),*(S+i),*(S+i+1),*(Ro+i),*(Ro+i+1),p,temp);
tP=tP*(*(temp+2));
}

zoep_P(*(P+m-2),*(P+m-1),*(S+m-2),*(S+m-1),*(Ro+m-2),*(Ro+m-1),p,temp);
tP=tP*(*temp);
for(i=m-1;i<2*m-3;i++)
{zoep_P(*(P+2*m-3-i),*(P+2*m-4-i),*(S+2*m-3-i),*(S+2*m-4-i),*(Ro+2*m-3-i),*(Ro+2*m-4-i),p,temp);
tP=tP*(*(temp+2));
}
}
*P_time=Ptime;
*P_para=tP;
*P_thita=pthita; //*
}

void S_1x_tpara(float *P,float *S,float *Ro,float *h,float x,int m,double *S_time,double *S_para,double *S_thita)//*
{double p,tS=1.0,*temp,Stime=0;
int i,it;
double sthita; //*
temp=calloc(4,sizeof(double));
p=S_p(h,P,S,x,m,0.000000001);
sthita=asin(p*(*S)); //*
for(it=0;it<m-1;it++)
{Stime=Stime+(*(h+it))/sqrt(1-p*p*(*(P+it))*(*(P+it)))/(*(P+it))+(*(h+it))/sqrt(1-p*p*(*(S+it))*(*(S+it)))/(*(S+it));
}
if(m==2)
{zoep_P(*P,*(P+1),*S,*(S+1),*Ro,*(Ro+1),p,temp);
tS=*(temp+1);
}
else
{ for(i=0;i<m-2;i++)
{zoep_P(*(P+i),*(P+i+1),*(S+i),*(S+i+1),*(Ro+i),*(Ro+i+1),p,temp);
tS=tS*(*(temp+2));
}

zoep_P(*(P+m-2),*(P+m-1),*(S+m-2),*(S+m-1),*(Ro+m-2),*(Ro+m-1),p,temp);
tS=tS*(*(temp+1));
for(i=m-1;i<2*m-3;i++)
{zoep_S(*(P+2*m-3-i),*(P+2*m-4-i),*(S+2*m-3-i),*(S+2*m-4-i),*(Ro+2*m-3-i),*(Ro+2*m-4-i),p,temp);
tS=tS*(*(temp+3));
}
}
*S_time=Stime;
*S_para=tS;
*S_thita=sthita; //*
}

void P_1x_allay_para(float *P,float *S,float *Ro,float *h,float x,int layn,double *P_time,double *P_para,double *P_thita) //*
{int m;
int i=0;
for(m=2;m<layn+1;m++)
{P_1x_tpara(P,S,Ro,h,x,m,P_time+i,P_para+i,P_thita+i); //*
i++;
}
}
void S_1x_allay_para(float *P,float *S,float *Ro,float *h,float x,int layn,double *S_time,double *S_para,double *S_thita)//*
{int m;
int i=0;
for(m=2;m<layn+1;m++)
{S_1x_tpara(P,S,Ro,h,x,m,S_time+i,S_para+i,S_thita+i);//*
i++;
}
}


main()
{double *P_time,*P_para,*P_thita,*S_time,*S_para,*S_thita;
float P[N]={1000,2000,2500,3000},S[N]={800,1800,2000,2000},Ro[N]={1,1,2,2.3},h[N]={1000,500,100,99999};
float x[20];
int nxo=20,i,j,pri;
for(j=0;j<20;j++)
{x[j]=500+j*50;}
P_time=calloc(10,sizeof(double));
P_para=calloc(10,sizeof(double));
P_thita=calloc(10,sizeof(double));//*
S_time=calloc(10,sizeof(double));
S_para=calloc(10,sizeof(double));
S_thita=calloc(10,sizeof(double));//*
for(i=0;i<nxo;i++)
{P_1x_allay_para(P,S,Ro,h,*(x+i),N,P_time,P_para,P_thita);//*
S_1x_allay_para(P,S,Ro,h,*(x+i),N,S_time,S_para,S_thita);//*
for(pri=0;pri<N-1;pri++)
{printf("P_time[%d]=%lf \n",pri,*(P_time+pri));
printf("P_para[%d]=%lf \n",pri,*(P_para+pri));
printf("P_thita[%d]=%lf \n",pri,*(P_thita+pri));
printf("S_time[%d]=%lf \n",pri,*(S_time+pri));
printf("S_para[%d]=%lf \n",pri,*(S_para+pri));
printf("S_thita[%d]=%lf \n",pri,*(S_thita+pri));
printf("***************\n");
}
printf("****************************************\n");
printf("****************************************\n");

}
}

[vfeeok1]

从网易邮箱添上去的网络图片好像没有显示。
那就把地址添上吧。
麻烦大家了。。
http://photo.163.com/wang.1989.jing/big/#aid=215935186&id=6638750377

[licaiyuren]

晕，你应该先自己分析好大概问题在哪边，你这贴代码这么多，谁看啊？？？？

[vfeeok1]

对，先说一下整个情况吧
我的目的是计算射线理论地震波的走时、界面的反射系数、地震波的出射的角度。
上面用*****************分隔的子函数，能单独实现一个功能。
第一个子函数是求射线参数的。
第二个子函数是求解一次佐普瑞兹方程。
第三个子函数是求解各参数的。
这三个子函数，我已经经过测试，确保没问题的。
***********************************************************************************

我估计问题是出现在main函数中。
情况是这样的，打印时屏幕上会出现部分计算结果，部分 fail。
但是在调试过程中，我所希望的值都能在alt+4的变量窗口显示。