运行环境Microsoft Visual C++ 6.0


#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));
    }
}