哪位好心人帮帮我，我的毕设！十万火急！！
因为初学c语言，很多地方不很明白，但现在面临毕业，很急，希望大家帮帮我。我想实现点光源的均匀照明，具体如下：
一点光源向四面八方发射光子，在其正下方有一个矩形区域，长度length，宽度width，我想设为4和3。原点在其中心（0,0,0，），设点光源在（0,0，h）处，在其上方有一个椭球面，x^2/a^2+y^2/b^2+(z-h)^2/c^=1,下面被平面z=z0截取，形成灯罩。发射大量光子，要求在矩形区域内光子均匀分布。
方案：
将矩形沿X轴方向以步长inter 等分成m 份，沿Y 轴方向以步长inter 等分成n 份，形成了m× n个等面积的长正方形网格，设在各个小格内的光子数构成数组gird[width/inter][length/inter]。
从光源p[3]=(0,0,h)射出一个光子,沿任意方向，方向向量i_l[3],求出与椭球的交点q[3]，判断是否射出椭球面，z坐标q[2]<z0,射出，反之，则求反射光线单位向量r_l[3],然后求新的交点，再判断。每射出一个光子，判断是否落在length*width的矩形范围内，若在范围内，统计总数total，并进一步判断落在那个小格内，统计各个小格内光子数gird[n][m],均匀照明即各个小格内光子数相等，理想情况各个小格内光子数ideal=total/3e4/4e4,(我将矩形分成3e4*4e4的网格)。
求目标函数objectfuction= ∑︳gird[][]-ideal︳，是由a,b,c,z0决定。
用模拟退火法，我首先假设a=0.05,b=0.04,c=0.04,z0=0.02,求出objectfuction，然后给a,b,c,z0加一个波动，再求objectfuction，如果objectfuction减小，新解被接受，否则以概率exp(-Δt′/T)接受作为新的当前解。 反复至最优解
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <time.h>
double a,b,c,z0;/*椭球曲面三个轴长及截面位置*/
float length=4.0,width=3.0;
float h=3.0;/*设定灯高*/
double p[3], i_l[3],n[3],q[3],r_l[3];/*p[3]椭圆内或椭圆上一点，i_l[3]直线的方向向量*/
/*求光线与椭球的交点*/
double intersect_point(p[3],i_l[3])
{
double a1,b1,c1;
double k,k1,k2,p;
a1 = i_l[0]*i_l[0]*b*b*c*c + i_l[1]*i_l[1]*a*a*c*c + i_l[2]*i_l[2]*a*a*b*b;
b1 = 2.0*i_l[0]*p[0]*b*b*c*c + 2.0*i_l[1]*p[1]*a*a*c*c + 2.0*i_l[2]*(p[2]-h)*a*a*b*b;
c1 = p[0]*p[0]*b*b*c*c + p[1]*p[1]*a*a*c*c + (p[2]-h)*(p[2]-h)*a*a*b*b- a*a*b*b*c*c; /*构造关于k的一元二次方程的参数a1k^2+b1k+c1=0*/
p = b1*b1 - 4*a1*c1;
k1 = (-b1 + sqrt(p))/a1 * 0.5;
k2 = (-b1 - sqrt(p))/a1 * 0.5;
k = (k1>k2)? k1:k2;
q[0] = k*i_l[0]+p[0];
q[1] = k*i_l[1]+p[1];
q[2] = k*i_l[2]+p[2];
return(q[3]);
}
/*求反射光线单位向量*/
double reject_light(p[3],i_l[3])
{
double t;
i_l[3]=i_l[3]/sqrt(i_l[0]*i_l[0]+i_l[1]*i_l[1]+i_l[2]*i_l[2]) ; /*入射光线单位向量*/
q[3]=intersect_point(p[3],i_l[3]);
n[3]={-2*q[0]/(a*a),-2*q[1]/(b*b),-2*(q[2]-h)/(c*c)}; /*椭球Q点出内侧法向量*/
n[3]=n[3]/sqrt(n[0]*n[0]+n[1]*n[1]+n[2]*n[2]); /*椭球Q点出内侧单位法向量*/
t =2*(i_l[0]*n[0]+i_l[1]*n[1]+i_l[2]*n[2]);
for(int i=0;i<3;i++)
r_l[i]=i_l[i]-t*n[i];
return(r_l[3]);
}
/*求各个方格内光子数*/
long count(a,b,c,z0)
{
double p[3]={0,0, h};/*光源位置*/
double x,y;
double inter=0.0001;
long m,n;
long gird[30000][40000];/*将目标矩形区域分成30000*40000的网格*/
long ideal;
long i,j,k=0;
time_t t;
for(i=1;i<=1e100;i++)/*发射10^100个光子*/
{
for(j=0;j<=2;j++) /*产生随机i_l[3]*/
{
srand(time(NULL));/*这里生成第一个种子值*/
i_l[i]=(double) rand()/RAND_MAX-0.5;
} /*随机产生一个数组,构成随机入射向量*/
q[3]=intersect_point( p[3],i_l[3]);
while(q[2]>z0)
{
k++;
if(k>10000)
break;
r_l[3]=reject_light(p[3],i_l[3]);
p[3]=q[3];
i_l[3]=r_l[3];
q[3]=intersect_point(p[3],i_l[3]);
}
if(k<=10000)
{
x=-p[2]*i_l[0]/i_l[2]+p[0];
y=-p[2]*i_l[1]/i_l[2]+p[1]; /*光子落在目标平面位置*/
if(-length/2<=x<=length/2&&-width/2<=y<=width/2)
{
m=floor((x+length/2)/inter);
n=floor((y+width/2)/inter);
gird[n][m]=gird[n][m]+1;/*求射在各方格内的光子数*/
}
}
}
return gird[30000][40000];
}
/*落在目标平面内的光子总数*/
long total(a,b,c,z0)
{
long num=0;
long i,j;
long gird[30000][40000]=count(a,b,c,z0);
for(i=0;i<3e4;i++)
{
for(j=0;j<4e4;j++)
num=num+gird[i][j];
}
return(num);
}
/*求光能利用率*/
float efficiency(a,b,c,z0)
{
float f=total(a,b,c,z0)/(1e100)*100%;
return(f);
}
/*求目标函数*/
double objectfunction(a,b,c,z0)
{
long sum=0;
long i,j;
long gird[30000[40000]=count(a,b,c,z0);
long ideal=floor(total(a,b,c,z0)/3e4/4e4);/*落在目标平面内光子实现均匀照明*/
for(i=0;i<3e4;i++)
{
for(j=0;j<4e4;j++)
{
gird[i][j]=abs(gird[i][j]-ideal);
sum=sum+gird[i][j];
}
}
return(sum);
}
/* 产生随机波动因子 */
double rnd()
{
double r;
r=(double) rand()/RAND_MAX;
return r;
}
/*模拟退火算法*/
main()
{
double objectfunction(a, b, c, z0);/*声明目标函数*/
long i;
int markovlength=1000000;/*马可夫链长度*/
double decayscale=0.95;/*衰减参数*/
double stepfactor=0.0002;/*Metropolis的步长*/
double temperature=100;/*初始温度*/
double tolerance=1e3;
double prea,preb,prec,prez0;
double nexta,nextb,nextc,nextz0;
double prebesta,prebestb,prebestc,prebestz0;
double besta,bestb,bestc,bestz0;
double acceptpoints=0.0;
time_t t;
prea=0.05;preb=0.04;prec=0.04;prez0=2.98;
srand((unsigned) time(&t));
prebesta=besta=prea;
prebestb=bestb=preb;
prebestc=bestc=prec;
prebestz0=bestz0=prez0;
do
{
temperature*=decayscale;
acceptpoints=0.0;
for(i=0;i<markovlength;i++)
{
do
{
nexta=prea+stepfactor*(rnd()-0.5);
nextb=preb+stepfactor*(rnd()-0.5);
nextc=prec+stepfactor*(rnd()-0.5);
nextz0=prez0+stepfactor*(rnd()-0.5);
}
while(efficiency(nexta,nextb,nextc,nextz0)<0.95);
if (objectfunction(besta,bestb,bestc,bestz0)>objectfunction(nexta,nextb,nextc,nextz0))
{
prebesta=besta;
prebestb=bestb;
prebestc=bestc;
prebestz0=bestz0;
besta=nexta;
bestb=nextb;
bestc=nextc;
bestz0=nextz0;
}
if (objectfunction(prea,preb,prec,prez0)-objectfunction(nexta,nextb,nextc,nextz0)>0)
{
prea=nexta;
preb=nextb;
prec=nextc;
prez0=nextz0;
acceptpoints++;
}
else
{
double change=-1*(objectfunction(nexta,nextb,nextc,nextz0)-objectfunction(prea,preb,prec,prez0))/temperature;
if(exp(change)>rnd())
{
prea=nexta;
preb=nextb;
prec=nextc;
prez0=nextz0;
acceptpoints++;
}
}
}
while(fabs(objectfunction(besta,bestb,bestc,bestz0)-objectfunction(prea,preb,prec,prez0))>tolerance);
printf("abs=%g\n",fabs(objectfunction(besta,bestb,bestc,bestz0)-objectfunction(prea,preb,prec,prez0)));
printf("%lf,%lf,%lf,%lf,%lf,%lf\n",prea,preb,prec,prez0,objectfunction(prea,preb,prec,prez0),temperature);
printf("The Min points is:%lf,%lf,%lf,%lf\n",besta,bestb,bestc,bestz0);
printf("The Min data is:%lf\n",objectfunction(besta,bestb,bestc,bestz0));
}
}
