StrayLight.cpp

// #include"StrayLight.h"
#include<iostream>
#include<string>
#include<fstream>
#include<cmath>
#include<string>
#include<ctime>
#include<iomanip>
#include<sstream>
#include<string.h>
#include"ephcom_wrapper.hpp"
using namespace std;

#define PI acos(-1)
#define Es 1367.0
#define Em 0.0044
#define Re 6371.0//地球半径
#define Re2 Re*Re
#define Swide 2.44//2023年4月1日结构设定为2.4425m，以后还可能变化，为方便计算设为2.44，不会有数量级误差
#define Slong 3.01//2023年4月1日结构设定为3.01m
#define Rd 1.1//2023年4月1日结构设定为1.128m，以后还可能变化，为方便计算设为1.1，不会有数量级误差
#define Plank 6.626e-34	//普朗克常数J*s
#define AU 149597870.7 //天文单位km
#define AU2 AU*AU

const double Es_fil[7] = {27.01477,86.82985,281.07176,255.6571,193.52574,200.77139,123.77403};
const double Em_fil[7] = {8.695317e-05,0.00027948159,0.0009046933,0.00082289045,0.00062290655,0.00064622833,0.00039839482};
const double PlanetsAbsMag[8][7]={2.231,0.782,-0.419,-0.773,-0.798,-0.718,-0.672,//水星
-3.367,-3.872,-4.410,-4.512,-4.398,-4.212,-4.117,//金星
1.130,0.101,-1.523,-2.520,-2.650,-2.503,-2.358,//火星
-8.528,-9.033,-9.402,-9.386,-9.047,-8.448,-8.421,//木星
-8.708,-9.213,-9.708,-9.764,-9.396,-8.608,-8.368,//土星
-6.346,-6.851,-7.197,-6.673,-5.336,-4.912,-5.317,//天王星
-5.933,-6.438,-6.900,-6.659,-5.470,-4.837,-5.081,//海王星
0.664,0.159,-0.719,-1.188,-1.196,-0.974,-0.891};//冥王星
const double EM0 = 2.8e-8;//0星等的照度值，公式：double E0 = 1365.08 / pow(100, 26.72 / 5);
//const double PlanetsAbsE[8][7]={3.58722e-09,1.36257e-08,4.11868e-08,5.70635e-08,5.83926e-08,5.42448e-08,5.19946e-08,
//6.2224e-07,9.90736e-07,1.62614e-06,1.78631e-06,1.60827e-06,1.35506e-06,1.24153e-06,
//9.88913e-09,2.55128e-08,1.13857e-07,2.85206e-07,3.21483e-07,2.80775e-07,2.45673e-07,
//7.21702e-05,0.00011491,0.00016142,0.000159059,0.000116401,6.70437e-05,6.5397e-05,
//8.5184e-05,0.000135631,0.000213973,0.000225298,0.000160531,7.76887e-05,6.22814e-05,
//9.67293e-06,1.54013e-05,2.11816e-05,1.30725e-05,3.81556e-06,2.58201e-06,3.74937e-06,
//6.61238e-06,1.05283e-05,1.61123e-05,1.2905e-05,4.31676e-06,2.40967e-06,3.01688e-06,
//1.519e-08,2.41857e-08,5.42948e-08,8.36292e-08,8.42477e-08,6.86686e-08,6.36148e-08};//水金火木土天海冥的绝对星等对应辐照度
const double PlanetsAbsE[8][7]={1.25183e-09,3.49052e-09,1.10171e-08,8.90266e-09,6.1621e-09,5.74084e-09,3.38622e-09,2.1714345e-07,2.5379912e-07,4.3497813e-07,2.7868908e-07,1.6971832e-07,1.4340872e-07,8.085651e-08,3.45101e-09,6.53567e-09,3.045562e-08,4.44959e-08,3.392569e-08,2.971499e-08,1.599982e-08,2.518527664e-05,2.94367667e-05,4.317848695e-05,2.481532382e-05,1.228368626e-05,7.0953786e-06,4.259072e-06,2.972670174e-05,3.474482318e-05,5.72357555e-05,3.514956436e-05,1.694059717e-05,8.22196406e-06,4.05615911e-06,3.37556695e-06,3.94539151e-06,5.66588062e-06,2.03948233e-06,4.0265094e-07,2.7325969e-07,2.4418275e-07,2.30752665e-06,2.69705688e-06,4.30990224e-06,2.0133531e-06,4.5554231e-07,2.5502078e-07,1.9647848e-07,5.30086e-09,6.1957e-09,1.452337e-08,1.304728e-08,8.89055e-09,7.26734e-09,4.143e-09};//水金火木土天海冥的绝对星等对应辐照度
const double PlanetsE_AU2[8][7]={28015446.582343537,78116167.85424984,246557514.80383125,199237356.53256097,137904849.78902856,128477299.08430137,75781955.27849847,4859566726.294486,5679903143.209061,9734602982.59355,6236928580.011572,3798214956.1848063,3209418842.794384,1809530078.5522048,77232047.71057059,146265173.2621918,681582206.1057857,995796902.2089956,759240802.5227964,665007277.5659081,358068244.0283927,563634475859.0736,658780794948.3726,966313938206.4478,555355108216.219,274903038143.9954,158791188059.1252,95315999534.81744,665269402821.2107,777572566648.6147,1280908901968.4758,786630481203.9672,379122482579.8927,184003633079.88455,90774905710.84904,75543577847.78088,88295979751.68683,126799705213.53505,45642641530.34039,9011135901.611265,6115421532.435761,5464693473.842928,51641345572.125656,60358846283.05587,96453556060.02104,45057881748.61057,10194819466.938034,5707243446.332199,4397094572.220029,118630824.75843115,138656722.37423745,325026105.3318717,291991877.7040579,198966307.20634192,162639637.65618783,92718373.92902023};//八个行星辐照度乘以一个天文单位的平方，便于事后计算实际场景辐照度
const double SunRatio[7] ={0.07209039,0.1050749,0.26816488,0.22462947,0.18818067,0.1100292,0.03183049};
const double MoonRatio[7]={0.07199216,0.10505516,0.26823942,0.22467527,0.18816909,0.11004074,0.03182816};
const int Planets[8]={MERCURY,VENUS,MARS,JUPITER,SATURN,URANUS,NEPTUNE,PLUTO};

extern "C"{
class XYZ//三维坐标类
{
public:
	double X;
	double Y;
	double Z;
	double Length;
	XYZ()
	{
		X = Y = Z = Length = 0;
	}
	XYZ(double x, double y, double z)
	{
		X = x, Y = y, Z = z, Length = Model();
	}
	XYZ(double a[3])
	{
	    X=a[0],Y=a[1],Z=a[2],Length=Model();
	}
	double Model()//模
	{
		return sqrt(X*X + Y*Y + Z*Z);
	}
	double Square()
	{
		return X*X + Y*Y + Z*Z;
	}
	double Cos(XYZ A)//夹角的cos值
	{
		return (X*A.X + Y*A.Y + Z*A.Z) / (Model()*A.Model());
	}
	double Angle(XYZ A)//夹角弧度
	{
		double a = X*A.X + Y*A.Y + Z*A.Z;
		double b = Model()*A.Model();
		double c = a / b;
		if (c > 1)return 0;
		else if (c < -1)return acos(-1);
		else return acos(c);
	}
	double Azimuth(XYZ A, XYZ py)//方位夹角弧度
	{
		double coefficient = (X*A.X + Y*A.Y + Z*A.Z) / (X*X + Y*Y + Z*Z);
		XYZ az = XYZ(A.X - X*coefficient, A.Y - Y*coefficient, A.Z - Z*coefficient);
		double c = az.Cos(py);
		return acos(c);
	}
	XYZ Cross(XYZ A)
	{
		return XYZ(Y*A.Z - Z*A.Y, Z*A.X - X*A.Z, X*A.Y - Y*A.X);
	}
	XYZ Normalize()
	{
		return XYZ(X / Length, Y / Length, Z / Length);
	}
	XYZ operator+(XYZ a)
	{
		return XYZ(X + a.X, Y + a.Y, Z + a.Z);
	}
	XYZ operator-(XYZ a)
	{
		return XYZ(X - a.X, Y - a.Y, Z - a.Z);
	}
	double operator*(XYZ a)
	{
		return X*a.X + Y*a.Y + Z*a.Z;
	}
};
class YZ//二维坐标类
{
public:
	double Y;
	double Z;
	YZ()
	{
		Y=Z=0;
	}
	YZ(double y,double z)
	{
		Y=y,Z=z;
	}
	double Model()//模
	{
		return sqrt(Y*Y+Z*Z);
	}
	double Square()
	{
		return Y*Y+Z*Z;
	}
	double Angle(YZ a)//夹角弧度
	{
		return acos((Y*a.Y+Z*a.Z)/(Model()*a.Model()));
	}
};

void Init(const char *deFn, const char *PSFFn, const char *RFn, const char *zolFn);
double PST(double Azimuth,double Angle);//73*45，行是方位角角5度一分，列是夹角2度一分，PST每个数占16字符，每行720+2字符

XYZ Spherical2Cartesian(double latitude,double longitude);//天球转三维，参数为弧度
void Cartesian2Spherical(XYZ car, double sph[2]);//所得sph是纬度经度的弧度
XYZ Ecliptic2Equatorial(XYZ Ecl);//黄道转赤道
XYZ Equatorial2Ecliptic(XYZ Equ);//赤道转黄道

double ShelterEarth(XYZ sat,XYZ star);//卫星是否会被地球遮挡
double ShelterPlate(XYZ sp,XYZ ob,XYZ py,double fai,double op,double sum);//亮星是否会被盖板遮挡

void PointSource(double jd, double sat[3], double ob[3], double py[3], double E[7], char* starMagFn);
void EarthShine(double ju, double sat[3], double ob[3], double py[3], double E[7]);
void Zodiacal(double ju, double ob[3], double E[7]);//琢磨半晌形参还是用赤道坐标系吧，内部再转换黄道！
void Calculate(double ju, double sat[3], double ob[3], double py[3], double E[7], char* starMagFn);

void ComposeY(double ob[3], double py1[3], double py2[3]);
void Watt2Photon(double E[7], double P[7]);//照度w/m2转光子数ph/(s*100um2)


double Pst[73][45];
const int Rrow = 180, Rcol = 288;
double R[Rrow][Rcol];
double ZL[37][19];


Ephcom ephcom = Ephcom();
void Init(const char *deFn, const char *PSFFn, const char *RFn, const char *zolFn)
{
//    string de405("DE405");
    string de405(deFn);
    ephcom.init(de405);
//    ifstream infile("PST");
    ifstream infile(PSFFn);
    string s;
    for (int i = 0; i < 73; i++)
    {
        getline(infile, s);
        stringstream ss(s);
        for (int j = 0; j < 45; j++)
        {
            ss >> Pst[i][j];
        }
    }
    infile.close();
//    infile.open("R");
    infile.open(RFn);
    for (int i = 0; i < Rrow; i++)
    {
        for (int j = 0; j < Rcol; j++)
        {
            infile >> R[i][j];
        }
    }
    infile.close();
//    infile.open("Zodiacal");
    infile.open(zolFn);
    for (int i = 0; i < 37; i++)
    {
        for (int j = 0; j < 19; j++)
        {
            infile >> ZL[i][j];
        }
    }
    infile.close();
}
double PST(double Azimuth,double Angle)//73*45，行是方位角角5度一分，列是夹角2度一分，PST每个数占16字符，每行720+2字符
{
    double daz=Azimuth*36/PI;
    double dan=Angle*90/PI;
    int az=daz;
    int an=dan;
    if(an==0)return 0;
    double pst=Pst[az][an-1]*(az+1-daz)*(an+1-dan)+Pst[az][an]*(az+1-daz)*(dan-an)+Pst[az+1][an-1]*(daz-az)*(an+1-dan)+Pst[az+1][an]*(daz-az)*(dan-an);
    return pst;
}

XYZ Spherical2Cartesian(double latitude,double longitude)//天球转三维，参数为弧度
{
    return XYZ(cos(latitude)*cos(longitude),cos(latitude)*sin(longitude),sin(latitude));
}
void Cartesian2Spherical(XYZ car, double sph[2])//所得sph是纬度经度的弧度
{
    sph[0] = asin(car.Z / car.Length);
    sph[1] = acos(car.X / sqrt(car.X*car.X + car.Y*car.Y));
    if (car.Y < 0)sph[1] = -sph[1];
}
XYZ Ecliptic2Equatorial(XYZ Ecl)//黄道转赤道
{
    double cosa = 0.917392;
    double sina =-0.397983;
    return XYZ(Ecl.X,Ecl.Y * cosa + Ecl.Z * sina,-Ecl.Y * sina + Ecl.Z * cosa);
}
XYZ Equatorial2Ecliptic(XYZ Equ)//赤道转黄道
{
    double cosa = 0.917392;
    double sina = 0.397983;
    return XYZ(Equ.X, Equ.Y * cosa + Equ.Z * sina, -Equ.Y * sina + Equ.Z * cosa);
}

double ShelterEarth(XYZ sat,XYZ star)//卫星是否会被地球遮挡
{
    double angle=star.Angle(sat);
    if(angle>PI/2)
    {
        double dis=sat.Model()*sin(angle);
        if(dis<Re)
        {
            return 0;
        }
        else
        {
            return 1;
        }
    }
    else
    {
        return 1;
    }
}
double ShelterPlate(XYZ sp,XYZ ob,XYZ py,double fai,double op,double sum)//亮星是否会被盖板遮挡
{
    double theta=sp.Angle(ob);
    double faisp=ob.Azimuth(sp,py);
    double faidiff=faisp-fai;
    if(theta<op||abs(faidiff)>=PI/2)
    {
        return 1;
    }
    else
    {
        double a=cos(theta);
        double b=sin(theta)*cos(faidiff);
        double c=sin(theta)*sin(faidiff);
        XYZ S=XYZ(a,b,c);
        double sumin=0;
        for(double y=-Rd; y<=Rd; y+=0.1)
        {
            for(double z=-sqrt(Rd*Rd-y*y); z<=sqrt(Rd*Rd-y*y); z+=0.1)
            {
                if((Slong*cos(op)-a/b*(Slong*sin(op)+Rd-y)>=0)&&(abs(z+(c*cos(op)*(Rd-y))/(b*cos(op)-a*sin(op)))<=Swide/2))
                {
                    sumin++;
                }
            }
        }
        return 1-sumin/sum;
    }
}

void PointSource(double ju, double sat[3], double ob[3], double py[3], double E[7], char* starMagFn)
{
    XYZ SAT = XYZ(sat[0], sat[1], sat[2]);
    XYZ OB = XYZ(ob[0], ob[1], ob[2]);
    XYZ PY = XYZ(py[0], py[1], py[2]);
    double pos[3], vel[3];
    ephcom.getPosVel(SUN, ju, pos, vel);
    XYZ sun = XYZ(pos[0], pos[1], pos[2]);
    double thetasun = OB.Angle(sun);
    if (thetasun < PI * 5 / 18)
    {
        cout << "观测方向与太阳夹角不可小于50度！目前只有" << thetasun * 180 / PI << "度！";
        return;
    }
    double faisun = OB.Azimuth(sun, PY);
    double fai, op;//盖板方位角，盖板打开角度-90°
    if (thetasun >= PI * 5 / 18 && thetasun < PI * 65 / 180)
    {
        op = 15 * PI / 180;
        fai = faisun;
    }
    else if (thetasun >= PI * 65 / 180 && thetasun < PI * 11 / 18)
    {
        op = 45 * PI / 180;
        fai = faisun;
    }
    else
    {
        op = 45 * PI / 180;
        fai = 0;
    }

    ephcom.getPosVel(MOON, ju, pos, vel);
    XYZ moon = XYZ(pos[0], pos[1], pos[2]);
    double thetamoon = OB.Angle(moon);
    double faimoon = OB.Azimuth(moon, PY);
    double moonphase = sun.Angle(moon) / PI;

    double E1[7] = { 0, 0, 0, 0, 0, 0, 0 };
    double pstsun, pstmoon;
    if (isfinite(faisun) == false){
        faisun = 0;
    }
    if (isfinite(thetasun) == false){
        thetasun = 0;
    }
    if (thetasun <= PI / 2)
    {
        pstsun = pow(10, PST(faisun, thetasun));
    }
    else
    {
        pstsun = 0;
    }
    if (isfinite(faimoon) == false){
        faimoon = 0;
    }
    if (isfinite(thetamoon) == false){
        thetamoon = 0;
    }
    if (thetamoon <= PI / 2 && thetamoon >= PI / 90)
    {
        pstmoon = pow(10, PST(faimoon, thetamoon));
    }
    else
    {
        pstmoon = 0;
    }
    double sum = 411;
    double Sheltersun, Sheltermoon;
    Sheltersun = ShelterEarth(SAT, sun)*ShelterPlate(sun, OB, PY, fai, op, sum);
    Sheltermoon = ShelterEarth(SAT, moon)*ShelterPlate(moon, OB, PY, fai, op, sum);
//    double Esun = Es*pstsun*Sheltersun, Emoon = Em*pstmoon*moonphase*Sheltermoon;
    for (int i = 0; i < 7; i++)
    {
        E1[i] = Es_fil[i]*pstsun*Sheltersun*SunRatio[i] + Em_fil[i]*pstmoon*moonphase*Sheltermoon*MoonRatio[i];
    }
    
    
    for (int i = 0; i < 2; i++)
    {
        ephcom.getPosVel(Planets[i], ju, pos, vel);
        XYZ planet = XYZ(pos[0], pos[1], pos[2]);
        double thetap = OB.Angle(planet);
        if (abs(thetap) < 0.008726646259971648){
            continue;
        }
        double faip = OB.Azimuth(planet, PY);
        if (isfinite(faip) == false){
            faip = 0;
        }
        if (isfinite(thetap) == false){
            thetap = 0;
        }
        double pstp = pow(10, PST(faip, thetap));
        double shelterp = ShelterEarth(SAT, planet)*ShelterPlate(planet, OB, PY, fai, op, sum);
        double planetphase = sun.Angle(planet) / PI;
        double proportion = pstp*shelterp*planetphase/planet.Square();
        for (int j = 0; j < 7; j++)
        {
            E1[j] += PlanetsE_AU2[i][j] * proportion;
        }
    }
    for (int i = 2; i < 8; i++)
    {
        ephcom.getPosVel(Planets[i], ju, pos, vel);
        XYZ planet = XYZ(pos[0], pos[1], pos[2]);
        double thetap = OB.Angle(planet);
        if (abs(thetap) < 0.008726646259971648){
            continue;
        }
        double faip = OB.Azimuth(planet, PY);
        if (isfinite(faip) == false){
            faip = 0;
        }
        if (isfinite(thetap) == false){
            thetap = 0;
        }
        double pstp = pow(10, PST(faip, thetap));
        double shelterp = ShelterEarth(SAT, planet)*ShelterPlate(planet, OB, PY, fai, op, sum);
        double proportion = pstp*shelterp/planet.Square();
        for (int j = 0; j < 7; j++)
        {
            E1[j] += PlanetsE_AU2[i][j] * proportion;
        }
    }

    string s;
    double Estar[7] = { 0, 0, 0, 0, 0, 0, 0 };
    ifstream infile(starMagFn);
    while (getline(infile, s))
    {
        stringstream ss(s);
        double x, y, z, d[7];
        ss >> x;
        ss >> y;
        ss >> z;
        for (int i = 0; i < 7; i++)
        {
            ss >> d[i];
        }
        XYZ star = XYZ(x, y, z);
        double theta = OB.Angle(star);
        if (abs(theta) < 0.008726646259971648){
            continue;
        }
        if (theta<PI / 2 && theta>PI / 90)
        {
            double faistar = OB.Azimuth(star, PY);
            if (isfinite(faistar) == false){
                faistar = 0;
            }
            if (isfinite(theta) == false){
                theta = 0;
            }
            double pst = pow(10, PST(faistar, theta));
            double shelter = ShelterEarth(SAT, star)*ShelterPlate(star,OB,PY,fai,op,sum);//太费时间了
            double proportion = pst*shelter;
            for (int i = 0; i < 7; i++)
            {
                Estar[i] += d[i] * proportion;
            }
        }
    }
    infile.close();

    // cout << "点源杂散光：";
    for (int i = 0; i < 7; i++)
    {
        E[i] = E1[i] + Estar[i];
        // cout << E[i] << " ";
    }
    // cout << endl;
}
void EarthShine(double ju, double sat[3], double ob[3], double py[3], double E[7])
{
    XYZ SAT = XYZ(sat[0], sat[1], sat[2]);
    XYZ OB = XYZ(ob[0], ob[1], ob[2]);
    XYZ PY = XYZ(py[0], py[1], py[2]);
    double pos[3], vel[3];
    ephcom.getPosVel(SUN, ju, pos, vel);
    XYZ sun = XYZ(pos[0], pos[1], pos[2]);
    ephcom.getPosVel(MOON, ju, pos, vel);
    XYZ moon = XYZ(pos[0], pos[1], pos[2]);
    double moonphase = sun.Angle(moon) / PI;
    double S[Rrow][Rcol];
    double wd = -PI / 2, jd = -PI, dwd = PI / 180, djd = PI / 150;
    double x[Rrow][Rcol], y[Rrow][Rcol], z[Rrow][Rcol], sunangle[Rrow][Rcol], moonangle[Rrow][Rcol];
    double Ese[Rrow][Rcol], Eme[Rrow][Rcol], Ee[Rrow][Rcol];
    for (int i = 0; i < Rrow; i++)
    {
        jd = -PI;
        for (int j = 0; j < Rcol; j++)
        {
            S[i][j] = Re2*cos(wd)*dwd*djd;
            x[i][j] = Re*cos(wd)*cos(jd);
            y[i][j] = Re*cos(wd)*sin(jd);
            z[i][j] = Re*sin(wd);
            XYZ v = XYZ(x[i][j], y[i][j], z[i][j]);
            sunangle[i][j] = v.Angle(sun);
            moonangle[i][j] = v.Angle(moon);
            if (sunangle[i][j] < PI / 2)
            {
                Ese[i][j] = Es*cos(sunangle[i][j])*R[i][j] * S[i][j];
                Eme[i][j] = 0;
            }
            else
            {
                Ese[i][j] = 0;
                if (moonangle[i][j] < PI / 2)
                {
                    Eme[i][j] = Em*moonphase*cos(moonangle[i][j])*R[i][j] * S[i][j];
                }
                else
                {
                    Eme[i][j] = 0;
                }
            }
            Ee[i][j] = Ese[i][j] + Eme[i][j];
            jd += djd;
        }
        wd += dwd;
    }

    double thetasun = OB.Angle(sun);
    double faisun = OB.Azimuth(sun, PY);
    double fai, op;
    if (thetasun < PI * 5 / 18)
    {
        cout << "观测方向与太阳夹角不可小于50度！目前只有" << thetasun * 180 / PI << "度！";
        return;
    }
    else if (thetasun >= PI * 5 / 18 && thetasun < PI * 65 / 180)
    {
        op = 15 * PI / 180;
        fai = faisun;
    }
    else if (thetasun >= PI * 65 / 180 && thetasun < PI * 11 / 18)
    {
        op = 45 * PI / 180;
        fai = faisun;
    }
    else
    {
        op = 45 * PI / 180;
        fai = 0;
    }

    OB = OB.Normalize();
    PY = PY.Normalize();
    XYZ zccst = OB.Cross(PY);
    XYZ sfa = XYZ(sin(op), -cos(op)*cos(fai), -cos(op)*sin(fai));//盖板法向量
    double sw = -1.1, sl = 0.1, dlw = 0.2, dlw2 = 0.04;//盖板z轴、y轴坐标和小块长度、小块面积
    int shieldW=Swide/dlw,shieldL=Slong/dlw;//盖板z轴\y轴方向，也就是宽\长方向，的小块数量；
    double xs[shieldW][shieldL], ys[shieldW][shieldL], zs[shieldW][shieldL], Eshield[shieldW][shieldL];//盖板小块坐标和接收照度
    double Rd2 = Rd*Rd;
    for (int w = 0; w < shieldW; w++)
    {
        for (int l = 0; l < shieldL; l++)
        {
            Eshield[w][l] = 0;
        }
    }
    for (int i = 0; i < Rrow; i++)
    {
        for (int j = 0; j < Rcol; j++)
        {
            XYZ v = XYZ(x[i][j] - SAT.X, y[i][j] - SAT.Y, z[i][j] - SAT.Z);//设施到地块
            if (XYZ(x[i][j], y[i][j], z[i][j])*v < 0)//判断地块是否为地球本身遮挡
            {
                XYZ vp = XYZ(v*OB, v*PY, v*zccst);//设施到地块转化为设施坐标系
                double Rv = v.Model();//设施到地块距离
                double Rv2 = v.Square();
                double coses = sfa*vp / Rv;//俩向量夹角的cos
                if (coses > 0)
                {
                    sw = -1.1;
                    for (int w = 0; w < shieldW; w++)
                    {
                        sl = 0.1;
                        for (int l = 0; l<shieldL; l++)
                        {
                            double y0 = Rd + sl*sin(op);
                            double z0 = sw;
                            double R0 = sqrt(y0*y0 + z0*z0);
                            double fai0 = asin(z0 / R0);
                            xs[w][l] = sl*cos(op);
                            ys[w][l] = R0*cos(fai0 + fai);
                            zs[w][l] = R0*sin(fai0 + fai);
                            double yjd = ys[w][l] - xs[w][l] * vp.Y / vp.X;
                            double zjd = zs[w][l] - xs[w][l] * vp.Z / vp.X;
                            if (yjd*yjd + zjd*zjd>Rd2)
                            {
                                Eshield[w][l] += Ee[i][j] * coses*dlw2 / 2 / PI / Rv2;
                            }
                            sl += dlw;
                        }
                        sw += dlw;
                    }
                }
            }
        }
    }

    double dlr = 0.1, dlt = 0.5*PI;//镜筒分块的半径和角度差
    int Mrow=Rd/dlr,Mcol=2*PI/dlt;
    double ye[Mrow][Mcol], ze[Mrow][Mcol];//镜筒分块的y坐标和z坐标，x必为0,故可省略
    double SM[Mrow][Mcol];//镜筒每个小块的面积
    double Econe[Mrow][Mcol];
    double Ealbedo1 = 0, Ealbedo2 = 0;//盖板反射进的地气光，直接射入的地气光
    for (int i = 0; i < Mrow; i++)
    {
        for (int j = 0; j < Mcol; j++)
        {
            ye[i][j] = dlr*(i+0.5)*cos(dlt*j);
            ze[i][j] = dlr*(i+0.5)*sin(dlt*j);
            SM[i][j]=dlr*(i+0.5)*dlt*dlr;//分块中心弧长乘以分块半径，约等于梯形面积
            Econe[i][j]=0;
            for (int w = 0; w < shieldW; w++)
            {
                for (int l = 0; l < shieldL; l++)
                {
                    XYZ ves = XYZ(xs[w][l], ys[w][l] - ye[i][j], zs[w][l] - ze[i][j]);//镜头分块到盖板小块
                    double tyj = acos(ves.X / ves.Length);
                    double fwj = acos(ves.Y / (sqrt(ves.Y*ves.Y + ves.Z*ves.Z)));
                    if (isfinite(tyj) == false){
                        tyj = 0;
                    }
                    if (isfinite(fwj) == false){
                        fwj = 0;
                    }
                    double ppp = pow(10, PST(fwj, tyj));
                    Econe[i][j] += Eshield[w][l] / 2 / PI / ves.Square()*0.03*pow(10, PST(fwj, tyj));
                }
            }
            Econe[i][j]*=SM[i][j];
            Ealbedo1+=Econe[i][j];
        }
    }
    Ealbedo1=Ealbedo1/Rd2/PI;

    for (int i = 0; i < Rrow; i++)
    {
        for (int j = 0; j < Rcol; j++)
        {
            XYZ v = XYZ(x[i][j] - SAT.X, y[i][j] - SAT.Y, z[i][j] - SAT.Z);//设施到地块
            double earthangle = OB.Angle(v);//观测方向与设施到地块方向的夹角
            double faie = OB.Azimuth(v, PY);//观测方向与设施到地块方向的方位角
            if (earthangle < PI / 2 && XYZ(x[i][j], y[i][j], z[i][j])*v < 0)
            {
                double faie2 = faie - fai;
                if (earthangle < op || abs(faie2)>PI / 2)
                {
                    Ealbedo2 += pow(10, PST(faie, earthangle))*Ee[i][j] / 2 / PI / v.Square();
                }
                else
                {
                    double D = 39961, J = 0;
                    for (double yn = -Rd; yn < Rd; yn += 0.01)
                    {
                        for (double zn = -Rd; zn < Rd; zn += 0.01)
                        {
                            if (yn*yn + zn*zn < Rd2)
                            {
                                double a = cos(earthangle);
                                double b = sin(earthangle)*cos(faie2);
                                double c = sin(earthangle)*sin(faie2);
                                if ((Slong*cos(op) >= a / b*(Slong*sin(op) - yn + Rd)) && (abs(zn + c*cos(op)*(zn - Rd) / (a*sin(op) - b*cos(op))) <= Swide / 2))
                                {
                                    J++;
                                }
                            }
                        }
                    }
                    Ealbedo2 += pow(10, PST(faie, earthangle))*Ee[i][j] / 2 / PI / v.Square()*(1 - J / D);
                }
            }
        }
    }
    double Ealbedo = Ealbedo1 + Ealbedo2;

    // cout << "地气杂散光：";
    for (int i = 0; i < 7; i++)
    {
        E[i] = Ealbedo*SunRatio[i];
        // cout << E[i] << " ";
    }
    // cout << endl;
}
void Zodiacal(double ju, double ob[3], double E[7])
{
    XYZ OB = XYZ(ob[0], ob[1], ob[2]);
    double pos[3], vel[3];
    ephcom.getPosVel(SUN, ju, pos, vel);
    XYZ sun = XYZ(pos[0], pos[1], pos[2]);
    double sph[2];
    OB = Equatorial2Ecliptic(OB);
    Cartesian2Spherical(OB, sph);
    double latitude = sph[0] / PI * 180, longitude = sph[1] / PI * 180;
    sun = Equatorial2Ecliptic(sun);
    Cartesian2Spherical(sun, sph);
    latitude = abs(latitude);
    longitude -= sph[1] / PI * 180;
    longitude = abs(longitude);
    if (longitude > 180)longitude = 360 - longitude;
    int lo, la;
    double lor, lar;
    lo = longitude / 5;
    lor = longitude / 5.0 - lo;
    la = latitude / 5;
    lar = latitude / 5.0 - la;
    double zl = ZL[lo][la] * (1 - lor)*(1 - lar) + ZL[lo][la + 1] * (1 - lor)*lar + ZL[lo + 1][la] * lor*(1 - lar) + ZL[lo + 1][la + 1] * lor*lar;
    zl*=PI / 784 * 1e-11;
    double u[7]= {132.3,106.944,110.541,43.8336,19.3065,6.54026,1.5624};
    // cout << "黄道弥散光：";
    for(int i=0; i<7; i++)
    {
        E[i]=zl*u[i];
        // cout << E[i] << " ";
    }
    // cout << endl;
}

void Calculate(double ju, double sat[3], double ob[3], double py[3], double E[7], char* starMagFn)
{
    double Ep[7],Ee[7],Ez[7];
    PointSource(ju, sat,ob,py,Ep, starMagFn);
    EarthShine(ju, sat,ob,py,Ee);
    Zodiacal(ju,ob,Ez);
    // cout << "总杂散光：";
    for(int i=0; i<7; i++)
    {
        E[i]=Ep[i]+Ee[i]+Ez[i];
        // cout<<E[i]<<" ";
    }
    // cout << endl;
}

void ComposeY(double ob[3], double py1[3], double py2[3])
{
    XYZ OB=XYZ(ob);
    XYZ Eclob = Equatorial2Ecliptic(OB);
    double sph[2];
    Cartesian2Spherical(Eclob, sph);
    if (abs(abs(sph[0]) - PI / 2) < 1e-6)
    {
        cout << "光轴不可指向极点";
        return;
    }
    double decy1 = sph[0] + PI / 2, decy2 = sph[0] - PI / 2;
    double ray1, ray2;
    if (decy1 > PI / 2)
    {
        decy1 = PI - decy1;
        ray1 = sph[1] + PI;
    }
    else ray1 = sph[1];
    if (decy2 > -PI / 2)ray2 = sph[1];
    else
    {
        decy2 = -decy2 - PI;
        ray2 = sph[1] + PI;
    }
    XYZ eclpy1 = Spherical2Cartesian(decy1, ray1);
    XYZ eclpy2 = Spherical2Cartesian(decy2, ray2);
    XYZ PY1 = Ecliptic2Equatorial(eclpy1);
    XYZ PY2 = Ecliptic2Equatorial(eclpy2);
    py1[0]=PY1.X,py1[1]=PY1.Y,py1[2]=PY1.Z;
    py2[0]=PY2.X,py2[1]=PY2.Y,py2[2]=PY2.Z;
}

void Watt2Photon(double E[7], double P[7])//照度w/m2转光子数ph/(s*100um2)
{
    double c = 3e8;
    double lamda[7]= {275e-9,354e-9,475e-9,622e-9,763e-9,900e-9,990e-9};
    double co=1e-10/Plank/c;//m2到100um2的系数
    //cout << "各谱段在单个像元的光子数：";
    for(int i=0; i<7; i++)
    {
        P[i]=E[i]*lamda[i]*co;
        //cout<<P[i]<<" ";
    }
    //cout << endl;
}

void Zodiacal1(double lonlat[2], double E[7])
{
	double latitude = abs(lonlat[1]);
	double longitude = abs(lonlat[0]);
	if (longitude > 180)longitude = 360 - longitude;
	int lo, la;
	double lor, lar;
	lo = longitude / 5;
	lor = longitude / 5.0 - lo;
	la = latitude / 5;
	lar = latitude / 5.0 - la;
	double zl = ZL[lo][la] * (1 - lor)*(1 - lar) + ZL[lo][la + 1] * (1 - lor)*lar + ZL[lo + 1][la] * lor*(1 - lar) + ZL[lo + 1][la + 1] * lor*lar;
	zl*=PI / 784 * 1e-11;
	double u[7]={132.3,106.944,110.541,43.8336,19.3065,6.54026,1.5624};
	// cout << "黄道弥散光：";
	for(int i=0;i<7;i++)
	{
		E[i]=zl*u[i];
		// cout << E[i] << " ";
	}
	// cout << endl;
}
}