Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

r.sun: Incorrect beam radiation for some steep inclined surfaces #2534

Merged
merged 13 commits into from
Mar 23, 2023
Merged
18 changes: 15 additions & 3 deletions raster/r.sun/main.c
Original file line number Diff line number Diff line change
Expand Up @@ -1711,6 +1711,7 @@ void calculate(double singleSlope, double singleAspect, double singleAlbedo,
double locTimeOffset;
double latitude, longitude;
double coslat;
bool shouldBeBestAM, isBestAM;


struct SunGeometryConstDay sunGeom;
Expand Down Expand Up @@ -1838,7 +1839,7 @@ void calculate(double singleSlope, double singleAspect, double singleAlbedo,
}
sunVarGeom.zmax = zmax;
shadowoffset_base = (j % (numRows)) * n * arrayNumInt;
#pragma omp parallel firstprivate(q1,tan_lam_l,z1,i,shadowoffset,longitTime,coslat,coslatsq,latitude,longitude,sin_phi_l,latid_l,sin_u,cos_u,sin_v,cos_v,lum,gridGeom,elevin,aspin,slopein,civiltime,linkein,albedo,latin,coefbh,coefdh,incidout,longin,horizon,beam_rad,insol_time,diff_rad,refl_rad,glob_rad,mapset,per,decimals,str_step)
#pragma omp parallel firstprivate(q1,tan_lam_l,z1,i,shadowoffset,longitTime,coslat,coslatsq,latitude,longitude,sin_phi_l,latid_l,sin_u,cos_u,sin_v,cos_v,lum,gridGeom,elevin,aspin,slopein,civiltime,linkein,albedo,latin,coefbh,coefdh,incidout,longin,horizon,beam_rad,insol_time,diff_rad,refl_rad,glob_rad,mapset,per,decimals,str_step,shouldBeBestAM,isBestAM)
{
#pragma omp for schedule(dynamic) \
firstprivate(sunGeom,sunVarGeom,sunSlopeGeom,sunRadVar) \
Expand Down Expand Up @@ -1952,8 +1953,19 @@ void calculate(double singleSlope, double singleAspect, double singleAlbedo,

q1 = gridGeom.sinlat * cos_u * sin_v +
gridGeom.coslat * sin_u;
tan_lam_l = -cos_u * cos_v / q1;
sunSlopeGeom.longit_l = atan(tan_lam_l);

if (q1 != 0.0) {
tan_lam_l = -cos_u * cos_v / q1;
sunSlopeGeom.longit_l = atan(tan_lam_l);
isBestAM = (tan_lam_l > 0);
} else {
sunSlopeGeom.longit_l = pihalf;
isBestAM = true;
}

shouldBeBestAM = (0.0 < sunSlopeGeom.aspect && sunSlopeGeom.aspect <= M_PI);
sunSlopeGeom.shift12hrs = (shouldBeBestAM != isBestAM);

sunSlopeGeom.lum_C31_l = cos(latid_l) * sunGeom.cosdecl;
sunSlopeGeom.lum_C33_l = sin_phi_l * sunGeom.sindecl;

Expand Down
10 changes: 7 additions & 3 deletions raster/r.sun/rsunlib.c
Original file line number Diff line number Diff line change
Expand Up @@ -341,6 +341,10 @@ double lumcline2(struct SunGeometryConstDay *sungeom,
/* func = cube; */
sunVarGeom->isShadow = 0;

double timeOffset = 0;
if (sunSlopeGeom->shift12hrs)
timeOffset = M_PI;

if (useShadow()) {
length = 0;

Expand Down Expand Up @@ -383,7 +387,7 @@ double lumcline2(struct SunGeometryConstDay *sungeom,
}
*/
s = sunSlopeGeom->lum_C31_l
* cos(-sungeom->timeAngle - sunSlopeGeom->longit_l)
* cos(-sungeom->timeAngle - sunSlopeGeom->longit_l + timeOffset)
+ sunSlopeGeom->lum_C33_l; /* Jenco */
}

Expand All @@ -409,7 +413,7 @@ double lumcline2(struct SunGeometryConstDay *sungeom,
}
*/
s = sunSlopeGeom->lum_C31_l
* cos(-sungeom->timeAngle - sunSlopeGeom->longit_l)
* cos(-sungeom->timeAngle - sunSlopeGeom->longit_l + timeOffset)
+ sunSlopeGeom->lum_C33_l; /* Jenco */
}
}
Expand All @@ -424,7 +428,7 @@ double lumcline2(struct SunGeometryConstDay *sungeom,
}
*/
s = sunSlopeGeom->lum_C31_l
* cos(-sungeom->timeAngle - sunSlopeGeom->longit_l)
* cos(-sungeom->timeAngle - sunSlopeGeom->longit_l + timeOffset)
+ sunSlopeGeom->lum_C33_l; /* Jenco */
}

Expand Down
2 changes: 1 addition & 1 deletion raster/r.sun/sunradstruct.h
Original file line number Diff line number Diff line change
Expand Up @@ -73,7 +73,7 @@ struct SunGeometryVarSlope
double lum_C33_l;
double slope;
double aspect;

bool shift12hrs;
};


Expand Down