trac https://trac.osgeo.org/grass/ticket/498: (merge from devbr6)

 lat= is no longer usable- chop it out;
 make latin= and longin= work, but they give us only ~1.3% speed boost);
 give some variables more understandable names


git-svn-id: https://svn.osgeo.org/grass/grass/trunk@38772 15284696-431f-4ddb-bdfa-cd5b030d7da7

raster/r.sun/main.c

@@ -126,11 +126,11 @@ long int shadowoffset = 0;
 int varCount_global = 0;
 int bitCount_global = 0;
 int arrayNumInt = 1;
-float **z = NULL, **o = NULL, **s = NULL, **li = NULL, **a = NULL, **la =
-    NULL, **longitArray, **cbhr = NULL, **cdhr = NULL;
+float **z = NULL, **o = NULL, **s = NULL, **li = NULL, **a = NULL,
+    **latitudeArray = NULL, **longitudeArray, **cbhr = NULL, **cdhr = NULL;
 double op, dp;
 double invstepx, invstepy;
-double sr_min = 24., sr_max = 0., ss_min = 24., ss_max = 0.;
+double sunrise_min = 24., sunrise_max = 0., sunset_min = 24., sunset_max = 0.;
 
 float **lumcl, **beam, **insol, **diff, **refl, **globrad;
 unsigned char *horizonarray = NULL;
@@ -141,10 +141,10 @@ double civilTime;
  */
 double xmin, xmax, ymin, ymax;
 double declin, step, dist;
-double li_max = 0., li_min = 100., al_max = 0., al_min = 1.0, la_max = -90.,
-    la_min = 90.;
+double linke_max = 0., linke_min = 100., albedo_max = 0., albedo_min = 1.0,
+    lat_max = -90., lat_min = 90.;
 double offsetx = 0.5, offsety = 0.5;
-char *tt, *lt;
+char *ttime;
 
 /*
  * double slope;
@@ -208,7 +208,7 @@ int main(int argc, char *argv[])
     struct
     {
 	struct Option *elevin, *aspin, *aspect, *slopein, *slope, *linkein,
-	    *lin, *albedo, *longin, *alb, *latin, *lat, *coefbh, *coefdh,
+	    *lin, *albedo, *longin, *alb, *latin, *coefbh, *coefdh,
 	    *incidout, *beam_rad, *insol_time, *diff_rad, *refl_rad,
 	    *glob_rad, *day, *step, *declin, *ltime, *dist, *horizon,
 	    *horizonstep, *numPartitions, *civilTime;
@@ -328,26 +328,16 @@ int main(int argc, char *argv[])
     parm.latin->required = NO;
     parm.latin->gisprompt = "old,cell,raster";
     parm.latin->description =
-	_("Name of the latitudes input raster map [decimal degrees]");
+	_("Name of input raster map containing latitudes [decimal degrees]");
     parm.latin->guisection = _("Input_options");
 
-    if (parm.latin->answer == NULL) {
-	parm.lat = G_define_option();
-	parm.lat->key = "lat";
-	parm.lat->type = TYPE_DOUBLE;
-	parm.lat->required = NO;
-	parm.lat->description =
-	    _("A single value of latitude [decimal degrees]");
-	parm.lat->guisection = _("Input_options");
-    }
-
     parm.longin = G_define_option();
     parm.longin->key = "longin";
     parm.longin->type = TYPE_STRING;
     parm.longin->required = NO;
     parm.longin->gisprompt = "old,cell,raster";
     parm.longin->description =
-	_("Name of the longitude input raster map [decimal degrees]");
+	_("Name of input raster map containing longitudes [decimal degrees]");
     parm.longin->guisection = _("Input_options");
 
     parm.coefbh = G_define_option();
@@ -553,20 +543,20 @@ int main(int argc, char *argv[])
     linkein = parm.linkein->answer;
     albedo = parm.albedo->answer;
     latin = parm.latin->answer;
+    longin = parm.longin->answer;
 
     civiltime = parm.civilTime->answer;
 
     if (civiltime != NULL) {
 	setUseCivilTime(TRUE);
-	longin = parm.longin->answer;
 
-	if (longin == NULL) {
-	    G_fatal_error(_("You must give the longitude raster if you use civil time"));
-
-	}
+	if (longin == NULL)
+	    G_fatal_error(
+		_("You must give the longitude raster if you use civil time"));
 
 	if (sscanf(parm.civilTime->answer, "%lf", &civilTime) != 1)
 	    G_fatal_error(_("Error reading civil time zone value"));
+
 	if (civilTime < -24. || civilTime > 24.)
 	    G_fatal_error(_("Invalid civil time zone value"));
 
@@ -614,7 +604,7 @@ int main(int argc, char *argv[])
 	G_fatal_error(_("If you use the horizon option you must also set the 'horizonstep' parameter."));
     }
 
-    tt = parm.ltime->answer;
+    ttime = parm.ltime->answer;
     if (parm.ltime->answer != NULL) {
 	if (insol_time != NULL)
 	    G_fatal_error(_("Time and insol_time are incompatible options"));
@@ -651,15 +641,6 @@ int main(int argc, char *argv[])
 	sscanf(parm.lin->answer, "%lf", &singleLinke);
     if (parm.albedo->answer == NULL)
 	sscanf(parm.alb->answer, "%lf", &singleAlbedo);
-    lt = parm.lat->answer;
-    /*
-     * if (parm.lat->answer != NULL)
-     * sscanf(parm.lat->answer, "%lf", &latitude);
-     */
-    /* HB 6/2008: why is the above commented out? instead of sscanf, maybe nicer to use:
-       G_scan_lat(parm.lat->answer, &latitude);
-       MN 2/2009: latitude doesn't exist! also G_scan_lat() does not exist in GRASS
-     */
 
     if (parm.slopein->answer == NULL)
 	sscanf(parm.slope->answer, "%lf", &singleSlope);
@@ -693,7 +674,9 @@ int main(int argc, char *argv[])
      * shadows pre-calculated, there is no problem. */
 
     if (saveMemory && useShadow() && (!useHorizonData()))
-	G_fatal_error(_("If you want to save memory and to use shadows, you must use pre-calculated horizons."));
+	G_fatal_error(
+	    _("If you want to save memory and to use shadows, "
+	      "you must use pre-calculated horizons."));
 
     if (parm.declin->answer == NULL)
 	declination = com_declin(day);
@@ -702,12 +685,7 @@ int main(int argc, char *argv[])
 	declination = -declin;
     }
 
-
-    /*
-     * if (lt != NULL)
-     * latitude = -latitude * deg2rad;
-     */
-    if (tt != 0) {
+    if (ttime != 0) {
 	/* Shadow for just one time during the day */
 	if (horizon == NULL) {
 	    arrayNumInt = 1;
@@ -723,7 +701,7 @@ int main(int argc, char *argv[])
 	}
     }
 
-    if (tt != NULL) {
+    if (ttime != NULL) {
 
 	tim = (timo - 12) * 15;
 	/* converting to degrees */
@@ -735,19 +713,22 @@ int main(int argc, char *argv[])
     }
 
     /* Set up parameters for projection to lat/long if necessary */
-
-    if (latin == NULL && lt == NULL && (G_projection() != PROJECTION_LL)) {
+    /*  we need to do this even if latin= and longin= were given because
+	rsunlib's com_par() needs iproj and oproj */
+    if ((G_projection() != PROJECTION_LL)) {
 	struct Key_Value *in_proj_info, *in_unit_info;
 
 	if ((in_proj_info = G_get_projinfo()) == NULL)
-	    G_fatal_error
-		(_("Can't get projection info of current location: please set latitude via 'lat' or 'latin' option!"));
+	    G_fatal_error(
+		_("Can't get projection info of current location"));
 
 	if ((in_unit_info = G_get_projunits()) == NULL)
-	    G_fatal_error(_("Can't get projection units of current location"));
+	    G_fatal_error(
+		_("Can't get projection units of current location"));
 
 	if (pj_get_kv(&iproj, in_proj_info, in_unit_info) < 0)
-	    G_fatal_error(_("Can't get projection key values of current location"));
+	    G_fatal_error(
+		_("Can't get projection key values of current location"));
 
 	G_free_key_value(in_proj_info);
 	G_free_key_value(in_unit_info);
@@ -760,8 +741,15 @@ int main(int argc, char *argv[])
 	    G_fatal_error(_("Unable to set up lat/long projection parameters"));
     }
 
+    if ((latin != NULL || longin != NULL) && (G_projection() == PROJECTION_LL))
+	G_warning(_("latin and longin raster maps have no effect when in a Lat/Lon location"));
+	/* true about longin= when civiltime is used? */
+    if ((latin == NULL && longin != NULL) || (latin != NULL && longin == NULL))
+	G_fatal_error(_("Both latin and longin raster maps must be given, or neither"));
+
 /**********end of parser - ******************************/
 
+
     if ((G_projection() == PROJECTION_LL))
 	ll_correction = TRUE;
 
@@ -785,8 +773,8 @@ int INPUT_part(int offset, double *zmax)
     FCELL *rast1 = NULL, *rast2 = NULL;
     static FCELL **horizonbuf;
     unsigned char *horizonpointer;
-    int fd1 = -1, fd2 = -1, fd3 = -1, fd4 = -1, fd5 = -1, fd6 = -1, fd7 =
-	-1, row, row_rev;
+    int fd1 = -1, fd2 = -1, fd3 = -1, fd4 = -1, fd5 = -1, fd6 = -1,
+	fd7 = -1, row, row_rev;
     static int *fd_shad;
     int fr1 = -1, fr2 = -1;
     int l, i, j;
@@ -877,10 +865,10 @@ int INPUT_part(int offset, double *zmax)
 
     if (latin != NULL) {
 	cell6 = Rast_allocate_f_buf();
-	if (la == NULL) {
-	    la = (float **)G_malloc(sizeof(float *) * (numRows));
+	if (latitudeArray == NULL) {
+	    latitudeArray = (float **)G_malloc(sizeof(float *) * (numRows));
 	    for (l = 0; l < numRows; l++)
-		la[l] = (float *)G_malloc(sizeof(float) * (n));
+		latitudeArray[l] = (float *)G_malloc(sizeof(float) * (n));
 	}
 	if ((mapset = G_find_raster2(latin, "")) == NULL)
 	    G_fatal_error(_("Raster map <%s> not found"), latin);
@@ -889,9 +877,9 @@ int INPUT_part(int offset, double *zmax)
 
     if (longin != NULL) {
 	cell7 = Rast_allocate_f_buf();
-	longitArray = (float **)G_malloc(sizeof(float *) * (numRows));
+	longitudeArray = (float **)G_malloc(sizeof(float *) * (numRows));
 	for (l = 0; l < numRows; l++)
-	    longitArray[l] = (float *)G_malloc(sizeof(float) * (n));
+	    longitudeArray[l] = (float *)G_malloc(sizeof(float) * (n));
 
 	if ((mapset = G_find_raster2(longin, "")) == NULL)
 	    G_fatal_error(_("Raster map <%s> not found"), longin);
@@ -932,7 +920,7 @@ int INPUT_part(int offset, double *zmax)
 	    fd_shad = (int *)G_malloc(sizeof(int) * arrayNumInt);
 	}
 	/*
-	 * if(tt != NULL)
+	 * if(ttime != NULL)
 	 * {
 	 * 
 	 * horizonbuf[0]=Rast_allocate_f_buf();
@@ -1041,16 +1029,16 @@ int INPUT_part(int offset, double *zmax)
 
 	    if (latin != NULL) {
 		if (!Rast_is_f_null_value(cell6 + j))
-		    la[rowrevoffset][j] = (float)cell6[j];
+		    latitudeArray[rowrevoffset][j] = (float)cell6[j];
 		else
-		    la[rowrevoffset][j] = UNDEFZ;
+		    latitudeArray[rowrevoffset][j] = UNDEFZ;
 	    }
 
 	    if (longin != NULL) {
 		if (!Rast_is_f_null_value(cell7 + j))
-		    longitArray[rowrevoffset][j] = (float)cell7[j];
+		    longitudeArray[rowrevoffset][j] = (float)cell7[j];
 		else
-		    longitArray[rowrevoffset][j] = UNDEFZ;
+		    longitudeArray[rowrevoffset][j] = UNDEFZ;
 	    }
 
 	    if (coefbh != NULL) {
@@ -1141,7 +1129,7 @@ int INPUT_part(int offset, double *zmax)
 		    z[i][j] = UNDEFZ;
 		if (albedo != NULL && a[i][j] == UNDEFZ)
 		    z[i][j] = UNDEFZ;
-		if (latin != NULL && la[i][j] == UNDEFZ)
+		if (latin != NULL && latitudeArray[i][j] == UNDEFZ)
 		    z[i][j] = UNDEFZ;
 		if (coefbh != NULL && cbhr[i][j] == UNDEFZ)
 		    z[i][j] = UNDEFZ;
@@ -1335,7 +1323,7 @@ void joules2(struct SunGeometryConstDay *sunGeom,
 
     com_par(sunGeom, sunVarGeom, gridGeom, latitude, longitude);
 
-    if (tt != NULL) {		/*irradiance */
+    if (ttime != NULL) {		/*irradiance */
 
 	s0 = lumcline2(sunGeom, sunVarGeom, sunSlopeGeom, gridGeom,
 		       horizonpointer);
@@ -1810,7 +1798,8 @@ void calculate(double singleSlope, double singleAspect, double singleAlbedo,
 	for (i = 0; i < n; i++) {
 
 	    if (useCivilTime()) {
-		longitTime = -longitArray[arrayOffset][i] / 15.;
+		/* sun travels 15deg per hour, so 1 TZ every 15 deg and 15 TZs * 24hrs = 360deg */
+		longitTime = -longitudeArray[arrayOffset][i] / 15.;
 	    }
 
 	    gridGeom.xg0 = gridGeom.xx0 = (double)i *gridGeom.stepx;
@@ -1841,42 +1830,49 @@ void calculate(double singleSlope, double singleAspect, double singleAlbedo,
 		}
 		if (linkein != NULL) {
 		    sunRadVar.linke = li[arrayOffset][i];
-		    li_max = AMAX1(li_max, sunRadVar.linke);
-		    li_min = AMIN1(li_min, sunRadVar.linke);
+		    linke_max = AMAX1(linke_max, sunRadVar.linke);
+		    linke_min = AMIN1(linke_min, sunRadVar.linke);
 		}
 		if (albedo != NULL) {
 		    sunRadVar.alb = a[arrayOffset][i];
-		    al_max = AMAX1(al_max, sunRadVar.alb);
-		    al_min = AMIN1(al_min, sunRadVar.alb);
+		    albedo_max = AMAX1(albedo_max, sunRadVar.alb);
+		    albedo_min = AMIN1(albedo_min, sunRadVar.alb);
 		}
 		if (latin != NULL) {
-		    latitude = la[arrayOffset][i];
-		    la_max = AMAX1(la_max, latitude);
-		    la_min = AMIN1(la_min, latitude);
+		    latitude = latitudeArray[arrayOffset][i];
+		    lat_max = AMAX1(lat_max, latitude);
+		    lat_min = AMIN1(lat_min, latitude);
 		    latitude *= deg2rad;
 		}
-		/* MN 2/2009: should it be?? 
-		   if (latin == NULL && lt == NULL && (G_projection() != PROJECTION_LL)) { 
-		 */
+		if (longin != NULL) {
+		    longitude = longitudeArray[arrayOffset][i];
+		    /* lon_max = AMAX1(lon_max, longitude); */
+		    /* lon_min = AMIN1(lon_min, longitude); */
+		    longitude *= deg2rad;
+		}
+
 		if ((G_projection() != PROJECTION_LL)) {
 
-		    longitude = gridGeom.xp;
-		    latitude = gridGeom.yp;
+		    if (latin == NULL || longin == NULL) {
+			/* if either is missing we have to calc both from current projection */
+			longitude = gridGeom.xp;
+			latitude = gridGeom.yp;
 
-		    if (pj_do_proj(&longitude, &latitude, &iproj, &oproj) < 0) {
-			G_fatal_error("Error in pj_do_proj");
-		    }
+			if (pj_do_proj(&longitude, &latitude, &iproj, &oproj) < 0) {
+			    G_fatal_error("Error in pj_do_proj");
+			}
 
-		    la_max = AMAX1(la_max, latitude);
-		    la_min = AMIN1(la_min, latitude);
-		    latitude *= deg2rad;
-		    longitude *= deg2rad;
+			lat_max = AMAX1(lat_max, latitude);
+			lat_min = AMIN1(lat_min, latitude);
+			latitude *= deg2rad;
+			longitude *= deg2rad;
+		    }
 		}
 		else {		/* ll projection */
 		    latitude = gridGeom.yp;
 		    longitude = gridGeom.xp;
-		    la_max = AMAX1(la_max, latitude);
-		    la_min = AMIN1(la_min, latitude);
+		    lat_max = AMAX1(lat_max, latitude);
+		    lat_min = AMIN1(lat_min, latitude);
 		    latitude *= deg2rad;
 		    longitude *= deg2rad;
 		}
@@ -1892,7 +1888,7 @@ void calculate(double singleSlope, double singleAspect, double singleAlbedo,
 		cos_v = cos(M_PI / 2 + sunSlopeGeom.aspect);
 		sin_v = sin(M_PI / 2 + sunSlopeGeom.aspect);
 
-		if (tt != NULL)
+		if (ttime != NULL)
 		    sunGeom.timeAngle = tim;
 
 		gridGeom.sinlat = sin(-latitude);
@@ -1912,10 +1908,10 @@ void calculate(double singleSlope, double singleAspect, double singleAlbedo,
 
 		if ((incidout != NULL) || someRadiation) {
 		    com_par_const(longitTime, &sunGeom, &gridGeom);
-		    sr_min = AMIN1(sr_min, sunGeom.sunrise_time);
-		    sr_max = AMAX1(sr_max, sunGeom.sunrise_time);
-		    ss_min = AMIN1(ss_min, sunGeom.sunset_time);
-		    ss_max = AMAX1(ss_max, sunGeom.sunset_time);
+		    sunrise_min = AMIN1(sunrise_min, sunGeom.sunrise_time);
+		    sunrise_max = AMAX1(sunrise_max, sunGeom.sunrise_time);
+		    sunset_min = AMIN1(sunset_min, sunGeom.sunset_time);
+		    sunset_max = AMAX1(sunset_max, sunGeom.sunset_time);
 
 		}
 
@@ -1985,7 +1981,7 @@ void calculate(double singleSlope, double singleAspect, double singleAlbedo,
 	    day);
     hist.edlinecnt = 2;
 
-    if (tt != NULL) {
+    if (ttime != NULL) {
 	sprintf(hist.edhist[hist.edlinecnt],
 		" Local (solar) time (decimal hr.):         %.4f", timo);
 	hist.edlinecnt++;
@@ -2000,17 +1996,12 @@ void calculate(double singleSlope, double singleAspect, double singleAlbedo,
 	    " Declination (rad):                        %f", -declination);
     hist.edlinecnt += 3;
 
-    if (lt != NULL)
-	sprintf(hist.edhist[hist.edlinecnt],
-		" Latitude (deg):                           %.4f",
-		-latitude * rad2deg);
-    else
-	sprintf(hist.edhist[hist.edlinecnt],
-		" Latitude min-max(deg):                    %.4f - %.4f",
-		la_min, la_max);
+    sprintf(hist.edhist[hist.edlinecnt],
+	    " Latitude min-max(deg):                    %.4f - %.4f",
+	    lat_min, lat_max);
     hist.edlinecnt++;
 
-    if (tt != NULL) {
+    if (ttime != NULL) {
 	sprintf(hist.edhist[hist.edlinecnt],
 		" Sunrise time (hr.):                       %.2f",
 		sunGeom.sunrise_time);
@@ -2024,16 +2015,16 @@ void calculate(double singleSlope, double singleAspect, double singleAlbedo,
     else {
 	sprintf(hist.edhist[hist.edlinecnt],
 		" Sunrise time min-max (hr.):               %.2f - %.2f",
-		sr_min, sr_max);
+		sunrise_min, sunrise_max);
 	sprintf(hist.edhist[hist.edlinecnt + 1],
 		" Sunset time min-max (hr.):                %.2f - %.2f",
-		ss_min, ss_max);
+		sunset_min, sunset_max);
 	sprintf(hist.edhist[hist.edlinecnt + 2],
 		" Time step (hr.):                          %.4f", step);
     }
     hist.edlinecnt += 3;
 
-    if (incidout != NULL || tt != NULL) {
+    if (incidout != NULL || ttime != NULL) {
 	sprintf(hist.edhist[hist.edlinecnt],
 		" Solar altitude (deg):                     %.4f",
 		sunVarGeom.solarAltitude * rad2deg);
@@ -2050,7 +2041,7 @@ void calculate(double singleSlope, double singleAspect, double singleAlbedo,
     else
 	sprintf(hist.edhist[hist.edlinecnt],
 		" Linke turbidity factor min-max:           %.1f-%.1f",
-		li_min, li_max);
+		linke_min, linke_max);
     hist.edlinecnt++;
 
     if (albedo == NULL)
@@ -2060,7 +2051,7 @@ void calculate(double singleSlope, double singleAspect, double singleAlbedo,
     else
 	sprintf(hist.edhist[hist.edlinecnt],
 		" Ground albedo min-max:                    %.3f-%.3f",
-		al_min, al_max);
+		albedo_min, albedo_max);
     hist.edlinecnt++;
 
     sprintf(hist.edhist[hist.edlinecnt],

raster/r.sun/r.sun.html

@@ -44,7 +44,8 @@ model accounts for sky obstruction by local relief features. Several solar
 parameters are saved in the resultant maps' history files, which may be viewed
 with the <a href="r.info.html">r.info</a> command.
 </p>
-<p>The solar incidence angle raster map <i>incidout</i> is computed specifying 
+<p>
+The solar incidence angle raster map <i>incidout</i> is computed specifying 
 elevation raster map <i>elevin</i>, aspect raster map <i>aspin</i>, slope 
 steepness raster map <i>slopin,</i> given the day <i>day</i> and local time
 <i>time</i>. There is no need to define latitude for locations with known
@@ -52,8 +53,7 @@ and defined projection/coordinate system (check it with the <a href="g.proj.html
 g.proj</a>
  command). If you have undefined projection, (x,y) system, etc. then the latitude
 can be defined explicitely for large areas by input raster map <i>latin</i>
- with interpolated latitude values or, for smaller areas, a single region 
-latitude value <i>lat</i> can be used instead. All input raster maps must
+ with interpolated latitude values. All input raster maps must
 be floating point (FCELL) raster maps. Null data in maps are excluded from
 the computation (and also speeding-up the computation), so each output raster
 map will contain null data in cells according to all input raster maps. The
@@ -63,7 +63,8 @@ The specified day <i>day</i> is the number of the day of the general year
 where January 1 is day no.1 and December 31 is 365. Time <i>time</i> must
 be a local (solar) time (i.e. NOT a zone time, e.g. GMT, CET) in decimal system,
 e.g. 7.5 (= 7h 30m A.M.), 16.1 = 4h 6m P.M.. </p>
-<p>Setting the solar declination <i>declin</i> by user is an option to override
+<p>
+Setting the solar declination <i>declin</i> by user is an option to override
 the value computed by the internal routine for the day of the year. The value
 of geographical latitude can be set as a constant for the whole computed
 region or, as an option, a grid representing spatially distributed values
@@ -71,15 +72,17 @@ over a large region. The geographical latitude must be also in decimal system
 with positive values for northern hemisphere and negative for southern one.
 In similar principle the Linke turbidity factor (<i>linkein</i>, <i>lin</i>
 ) and ground albedo (<i>albedo</i>, <i>alb</i>) can be set. </p>
-<p>Besides clear-sky radiations, the user can compute a real-sky radiation (beam,
+<p>
+Besides clear-sky radiations, the user can compute a real-sky radiation (beam,
 diffuse) using <i>coefbh</i> and <i>coefdh </i>input raster maps defining
 the fraction of the respective clear-sky radiations reduced by atmospheric
 factors (e.g. cloudiness). The value is between 0-1. Usually these
 coefficients can be obtained from a long-terms meteorological measurements.</p>
-<p>The solar irradiation or irradiance raster maps <i>beam_rad</i>, <i>diff_rad</i>
-, <i>refl_rad</i> are computed for a given day <i>day,</i> latitude <i>lat
-(latin), </i>elevation <i>elevin</i>, slope <i>slopein</i> and aspect <i>
-aspin</i> raster maps. For convenience, the output raster given as <i>glob_rad</i>
+<p>
+The solar irradiation or irradiance raster maps <i>beam_rad</i>, <i>diff_rad</i>,
+<i>refl_rad</i> are computed for a given day <i>day,</i> latitude <i>latin</i>,
+elevation <i>elevin</i>, slope <i>slopein</i> and aspect <i>aspin</i> raster maps.
+For convenience, the output raster given as <i>glob_rad</i>
 will output the sum of the three radiation components. The program uses 
 the Linke atmosphere turbidity factor and ground albedo coefficient. 
 A default, single value of Linke factor is <i>lin</i>=3.0 and 

raster/r.sun/rsunlib.c

@@ -183,7 +183,6 @@ void com_par(struct SunGeometryConstDay *sungeom,
     costimeAngle = cos(sungeom->timeAngle);
 
 
-
     lum_Lx = -sungeom->lum_C22 * sin(sungeom->timeAngle);
     lum_Ly = sungeom->lum_C11 * costimeAngle + sungeom->lum_C13;
     sunVarGeom->sinSolarAltitude =
@@ -240,7 +239,7 @@ void com_par(struct SunGeometryConstDay *sungeom,
     inputAngle = (inputAngle >= pi2) ? inputAngle - pi2 : inputAngle;
 
 
-    delt_lat = -0.0001 * cos(inputAngle);	/* Arbitrary small distance in latitude */
+    delt_lat = -0.0001 * cos(inputAngle);  /* Arbitrary small distance in latitude */
     delt_lon = 0.0001 * sin(inputAngle) / cos(latitude);
 
     newLatitude = (latitude + delt_lat) * rad2deg;
@@ -259,7 +258,6 @@ void com_par(struct SunGeometryConstDay *sungeom,
     delt_dist = sqrt(delt_east * delt_east + delt_nor * delt_nor);
 
 
-
     sunVarGeom->stepsinangle = gridGeom->stepxy * delt_nor / delt_dist;
     sunVarGeom->stepcosangle = gridGeom->stepxy * delt_east / delt_dist;