// #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();
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]);
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]);
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];
string de405("DE405");
Ephcom ephcom(de405);
void Init()
{
ifstream infile("PST");
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");
for (int i = 0; i < Rrow; i++)
{
for (int j = 0; j < Rcol; j++)
{
infile >> R[i][j];
}
}
infile.close();
infile.open("Zodiacal");
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])
{
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 (thetasun <= PI / 2)
{
pstsun = pow(10, PST(faisun, thetasun));
}
else
{
pstsun = 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);
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);
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("BrightGaia_with_csst_mag");
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);
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)));
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])
{
double Ep[7],Ee[7],Ez[7];
PointSource(ju, sat,ob,py,Ep);
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;
}
}