set term png set out 'r554_170.png' set grid set size 1.8 set xtics 10 set ytics 0.1 # A and B should be _k, i.e. wavelenght-dependent: A0 = -1.889 A1 = -1.628 A2 = 0.438 A3 = -0.235 #553.8 A0_554=-2.12504 A1_554=-1.65970 A2_554=0.38409 A3_554=-0.20655 #486.9 A0_487=-2.23351 A1_487=-1.68573 A2_487= 0.37632 A3_487=-0.19877 # 665.1 A0_665=-1.88914 A1_665=-1.58096 A2_665=0.30477 A3_665=-0.17908 B1_487=0.03881 B2_487=0.01566 B3_487=-0.00555 B1_665=0.04415 B2_665=0.00983 B3_665=-0.00389 B1 = 0.0425 B2 = 0.0132 B3 = -0.005 #553.8 B1_554= 0.04052 B2_554= 0.01009 B3_554= -0.00388 #C1 = 0.0003 #C2 = -0.0014 #C3 = 0.0010 #C4 = 0.0006 C1 = 0.00034115 # 0.0003 C2 = -0.0013425 # -0.0014 C3 = 0.00095906 # 0.0010 C4 = 0.00066229 # 0.0006 # D too, _k, i.e. wavelenght-dependent: D1 = 0.389 D2 = -0.148 D3 = -0.0035 D1_554= 0.37206 D2_554= -0.10745 D3_554= 0.00347 D1_487=0.36590 D2_487=-0.08945 D3_487=0.00678 D1_665=0.37141 D2_665=-0.13514 D3_665=0.01248 #P1 = 3.08 #P2 = 12.10 #P3 =-43.48 #P4 = 18.73 P1 = 4.06054 # 3.08 P2 = 12.8802 # 12.10 P3 =-30.5858 # -43.48 P4 = 16.7498 # 18.73 # t and f are librations (terr. lat., long.) / 1 degree # F is selen. long. of the Sun / 1 rad # g is the absolute value of phase angle / 1 rad g(x) = abs(pi*G(x)/180) G(x) = abs(sqrt(x**2 + t**2)) F(x) = pi*x/180 t = 0 f = 0 lnA_554(x) = \ A0_554 + A1_554*g(x) + A2_554*(g(x)**2) + A3_554*(g(x)**3) + \ B1_554*F(x) + B2_554*(F(x)**3) + B3_554*(F(x)**5) + \ C1*t + C2*f + C3*F(x)*t + C4*F(x)*f + \ D1_554*exp(-G(x)/P1) + D2_554*exp(-G(x)/P2) + D3_554*cos((G(x)-P3)/P4) logA_554(x) = lnA_554(x) / log(10) #+ 0.11 lnA_487(x) = \ A0_487 + A1_487*g(x) + A2_487*(g(x)**2) + A3_487*(g(x)**3) + \ B1_487*F(x) + B2_487*(F(x)**3) + B3_487*(F(x)**5) + \ C1*t + C2*f + C3*F(x)*t + C4*F(x)*f + \ D1_487*exp(-G(x)/P1) + D2_487*exp(-G(x)/P2) + D3_487*cos((G(x)-P3)/P4) logA_487(x) = lnA_487(x) / log(10) + 0.05 lnA_665(x) = \ A0_665 + A1_665*g(x) + A2_665*(g(x)**2) + A3_665*(g(x)**3) + \ B1_665*F(x) + B2_665*(F(x)**3) + B3_665*(F(x)**5) + \ C1*t + C2*f + C3*F(x)*t + C4*F(x)*f + \ D1_665*exp(-G(x)/P1) + D2_665*exp(-G(x)/P2) + D3_665*cos((G(x)-P3)/P4) logA_665(x) = lnA_665(x) / log(10) + -0.11 lnA(x) = \ A0 + A1*g(x) + A2*(g(x)**2) + A3*(g(x)**3) + \ B1*F(x) + B2*(F(x)**3) + B3*(F(x)**5) + \ C1*t + C2*f + C3*F(x)*t + C4*F(x)*f + \ D1*exp(-G(x)/P1) + D2*exp(-G(x)/P2) + D3*cos((G(x)-P3)/P4) logA(x) = lnA(x) / log(10) # (this would hold for t=5:) lnA5(x) = \ A0 + A1*g(x) + A2*(g(x)**2) + A3*(g(x)**3) + \ B1*F(x) + B2*(F(x)**3) + B3*(F(x)**5) + \ C1*(t+5) + C2*f + C3*F(x)*(t+5) + C4*F(x)*f + \ D1*exp(-G(x)/P1) + D2*exp(-G(x)/P2) + D3*cos((G(x)-P3)/P4) logA5(x) = lnA5(x) / log(10) fai(x) = 0.026 * abs(x) + 4.0E-9 * (x**4) logAold(x) = -0.95 - fai(x) / 2.5 set title 'log10 of lunar lambertian reflectance ( pi x flux density from / to Moon)' set xlabel 'selenographic long. of the Sun / 1 degree (positive for waxing Moon)' plot [-170:170] \ logA_487(x) title \ '0.05 + 486.9 nm ROLO formula, valid in (-90,90), here for zero libration' \ with lines 3, \ logA_554(x) title \ '553.8 nm ROLO formula, valid in (-90,90), here for zero libration' \ with lines 2, \ logA_665(x) title \ '-0.11 + 665.1 nm ROLO formula, valid in (-90,90), here for zero libration' \ with lines 1, \ logAold(x) title 'Old formula: -0.95 -(0.026 * abs(x) + 4.0E-9 * (x**4))/2.5' \ with lines 0 # logA(-9090 ? x : 1/0) title '' with lines 0, \ # Moon: +0.23 + 5*log10(r*R) + 0.026 * FV + 4.0E-9 * FV**4 # acc. to # Linkname: Computing planetary positions # URL: http://hotel04.ausys.se/pausch/comp/ppcomp.html#15