Merge pull request #60 from MennoVeerman/main

Fix bugs in tilting aerosols concentrations in tica mode.

Author: MennoVeerman

include_rt_kernels/raytracer_kernels_bw.h

@@ -27,7 +27,7 @@ constexpr int bw_kernel_grid = 256;
 // sun has a half angle of .266 degrees
 constexpr Float cos_half_angle = Float(0.9999891776066407); // cos(half_angle);
 constexpr Float sun_solid_angle = Float(6.799910294339209e-05); // 2.*M_PI*(1-cos_half_angle);
-constexpr Float sun_solid_angle_reciprocal = Float(14706.07635563193); 
+constexpr Float sun_solid_angle_reciprocal = Float(14706.07635563193);
 
 
 struct Grid_knull
@@ -40,7 +40,7 @@ struct Grid_knull
 struct Camera
 {
     Vector<Float> position;
-    bool fisheye = true;
+    int cam_type; // (0: fisheye, 1: rectangular, 2: top-of-atmosphere radiances)
 
     // rotation matrix for fisheye version - we need to do implement this in a nice way at some point
     Vector<Float> mx;
@@ -66,17 +66,17 @@ struct Camera
 
     void setup_normal_camera(const Camera camera)
     {
-        if (!fisheye)
+        if (camera.cam_type != 0)
         {
             const Vector<Float> dir_tmp = {1, 0, 0};
+            const Vector<Float> dir_up = {0, 0, 1};
+
             const Vector<Float> dir_cam = normalize(Vector<Float>({dot(camera.mx,dir_tmp), dot(camera.my,dir_tmp), dot(camera.mz,dir_tmp)}));
-            const Vector<Float> dir_up_tmp = {0, 0, 1};
-            const Vector<Float> dir_up = normalize(Vector<Float>({dot(camera.mx,dir_up_tmp), dot(camera.my,dir_up_tmp), dot(camera.mz,dir_up_tmp)}));
 
+            cam_height = normalize(Vector<Float>({dot(camera.mx, dir_up), dot(camera.my,dir_up), dot(camera.mz,dir_up)}));
             cam_width = Float(-1) * normalize(cross(dir_cam, dir_up));
-            cam_height = normalize(cross(dir_cam, cam_width));
             cam_depth = dir_cam / tan(fov/Float(180)*M_PI/Float(2.));
-        
+
             if (camera.nx > camera.ny)
                 cam_height = cam_height * Float(camera.ny)/Float(camera.nx);
             else if (camera.ny > camera.nx)

include_tilt/tilt_utils.h

@@ -45,35 +45,35 @@ void post_process_output(const std::vector<ColumnResult>& col_results,
         const bool switch_liq_cloud_optics,
         const bool switch_ice_cloud_optics);
 
-void restore_bkg_profile(const int n_x, const int n_y, 
+void restore_bkg_profile(const int n_x, const int n_y,
                       const int n_full,
-                      const int n_tilt, 
-                      const int bkg_start, 
+                      const int n_tilt,
+                      const int bkg_start,
                       std::vector<Float>& var,
                       std::vector<Float>& var_w_bkg);
 
-void restore_bkg_profile_bundle(const int n_col_x, const int n_col_y, 
-    const int n_lay, const int n_lev, 
-    const int n_lay_tot, const int n_lev_tot, 
+void restore_bkg_profile_bundle(const int n_col_x, const int n_col_y,
+    const int n_lay, const int n_lev,
+    const int n_lay_tot, const int n_lev_tot,
     const int n_z_in, const int n_zh_in,
-    const int bkg_start_z, const int bkg_start_zh, 
-    Array<Float,2>* p_lay_copy, Array<Float,2>* t_lay_copy, Array<Float,2>* p_lev_copy, Array<Float,2>* t_lev_copy, 
-    Array<Float,2>* lwp_copy, Array<Float,2>* iwp_copy, Array<Float,2>* rel_copy, Array<Float,2>* dei_copy, Array<Float,2>* rh_copy, 
+    const int bkg_start_z, const int bkg_start_zh,
+    Array<Float,2>* p_lay_copy, Array<Float,2>* t_lay_copy, Array<Float,2>* p_lev_copy, Array<Float,2>* t_lev_copy,
+    Array<Float,2>* lwp_copy, Array<Float,2>* iwp_copy, Array<Float,2>* rel_copy, Array<Float,2>* dei_copy, Array<Float,2>* rh_copy,
     Gas_concs& gas_concs_copy, Aerosol_concs& aerosol_concs_copy,
-    Array<Float,2>* p_lay, Array<Float,2>* t_lay, Array<Float,2>* p_lev, Array<Float,2>* t_lev, 
-    Array<Float,2>* lwp, Array<Float,2>* iwp, Array<Float,2>* rel, Array<Float,2>* dei, Array<Float,2>* rh, 
-    Gas_concs& gas_concs, Aerosol_concs& aerosol_concs, 
+    Array<Float,2>* p_lay, Array<Float,2>* t_lay, Array<Float,2>* p_lev, Array<Float,2>* t_lev,
+    Array<Float,2>* lwp, Array<Float,2>* iwp, Array<Float,2>* rel, Array<Float,2>* dei, Array<Float,2>* rh,
+    Gas_concs& gas_concs, Aerosol_concs& aerosol_concs,
     std::vector<std::string> gas_names, std::vector<std::string> aerosol_names,
     bool switch_liq_cloud_optics, bool switch_ice_cloud_optics, bool switch_aerosol_optics
     );
 
-void compress_columns_weighted_avg(const int n_x, const int n_y,  
-                      const int n_out, 
+void compress_columns_weighted_avg(const int n_x, const int n_y,
+                      const int n_out,
                       const int n_tilt,
                       const int compress_lay_start_idx,
                       std::vector<Float>& var, std::vector<Float>& var_weighting);
 
-void compress_columns_p_or_t(const int n_x, const int n_y, 
+void compress_columns_p_or_t(const int n_x, const int n_y,
                       const int n_out_lay,  const int n_tilt,
                       const int compress_lay_start_idx,
                       std::vector<Float>& var_lev, std::vector<Float>& var_lay);
@@ -82,16 +82,16 @@ void tilt_fields(const int n_z_in, const int n_zh_in, const int n_col_x, const i
     const int n_z_tilt, const int n_zh_tilt, const int n_col,
     const Array<Float,1> zh, const Array<Float,1> z,
     const Array<Float,1> zh_tilt, const Array<ijk,1> path,
-    Array<Float,2>* p_lay_copy, Array<Float,2>* t_lay_copy, Array<Float,2>* p_lev_copy, Array<Float,2>* t_lev_copy, 
-    Array<Float,2>* rh_copy, 
+    Array<Float,2>* p_lay_copy, Array<Float,2>* t_lay_copy, Array<Float,2>* p_lev_copy, Array<Float,2>* t_lev_copy,
+    Array<Float,2>* rh_copy,
     Gas_concs& gas_concs_copy, const std::vector<std::string> gas_names,
     Aerosol_concs& aerosol_concs_copy, const std::vector<std::string> aerosol_names, const bool switch_aerosol_optics
     );
 
 void compress_fields(const int compress_lay_start_idx, const int n_col_x, const int n_col_y,
     const int n_z_in, const int n_zh_in,  const int n_z_tilt,
-    Array<Float,2>* p_lay_copy, Array<Float,2>* t_lay_copy, Array<Float,2>* p_lev_copy, Array<Float,2>* t_lev_copy, 
-    Array<Float,2>* rh_copy, 
+    Array<Float,2>* p_lay_copy, Array<Float,2>* t_lay_copy, Array<Float,2>* p_lev_copy, Array<Float,2>* t_lev_copy,
+    Array<Float,2>* rh_copy,
     Gas_concs& gas_concs_copy, std::vector<std::string> gas_names,
     Aerosol_concs& aerosol_concs_copy, std::vector<std::string> aerosol_names, const bool switch_aerosol_optics);
 
@@ -114,11 +114,11 @@ void tica_tilt(const Float sza, const Float azi,
     const int n_col_x, const int n_col_y, const int n_col,
     const int n_lay, const int n_lev, const int n_z_in, const int n_zh_in ,
     Array<Float,1> xh, Array<Float,1> yh, Array<Float,1> zh, Array<Float,1> z,
-    Array<Float,2> p_lay, Array<Float,2> t_lay, Array<Float,2> p_lev, Array<Float,2> t_lev, 
-    Array<Float,2> lwp, Array<Float,2> iwp, Array<Float,2> rel, Array<Float,2> dei, Array<Float,2> rh, 
+    Array<Float,2> p_lay, Array<Float,2> t_lay, Array<Float,2> p_lev, Array<Float,2> t_lev,
+    Array<Float,2> lwp, Array<Float,2> iwp, Array<Float,2> rel, Array<Float,2> dei, Array<Float,2> rh,
     Gas_concs gas_concs, Aerosol_concs aerosol_concs,
-    Array<Float,2>& p_lay_out, Array<Float,2>& t_lay_out, Array<Float,2>& p_lev_out, Array<Float,2>& t_lev_out, 
-    Array<Float,2>& lwp_out, Array<Float,2>& iwp_out, Array<Float,2>& rel_out, Array<Float,2>& dei_out, Array<Float,2>& rh_out, 
-    Gas_concs& gas_concs_out, Aerosol_concs aerosol_concs_out,
+    Array<Float,2>& p_lay_out, Array<Float,2>& t_lay_out, Array<Float,2>& p_lev_out, Array<Float,2>& t_lev_out,
+    Array<Float,2>& lwp_out, Array<Float,2>& iwp_out, Array<Float,2>& rel_out, Array<Float,2>& dei_out, Array<Float,2>& rh_out,
+    Gas_concs& gas_concs_out, Aerosol_concs& aerosol_concs_out,
     std::vector<std::string> gas_names, std::vector<std::string> aerosol_names,
     bool switch_cloud_optics, bool switch_liq_cloud_optics, bool switch_ice_cloud_optics, bool switch_aerosol_optics);

python/image_from_xyz.py

@@ -72,8 +72,8 @@
 
     # reshape RGB array
     RGB = RGB.reshape((3, ny, nx))
-    RGB = RGB.swapaxes(0,2).swapaxes(0,1)[::-1,:]
-    
+    RGB = RGB.swapaxes(0,2).swapaxes(0,1)[:,:]
+
     # plot
     fig,ax = pl.subplots(figsize=(nx/args.dpi, ny/args.dpi), frameon=False)
     ax.grid(False)

python/set_virtual_camera.py

@@ -17,8 +17,7 @@
     "yaw": "Horizontal direction of camera, 0: east, 90: north, 180/-180: weast, -90/270: south",
     "pitch": "Vertical direction of camera, -90: vertical upward, 0: horizontal, 90: vertically downward.",
     "roll": "Roll of camera over the direction of camera",
-    "fisheye": "Whether to use fisheye lens (1) or a rectangular/square lens (0)",
-    "f_zoom": "Zoom factor (if fisheye=1)",
+    "cam_type": "camera type (0: fisheye lens, 1: rectangular lens, 2: top-of-atmosphere upwelling radiance)",
     "fov": "Field of view (if fisheye=0)",
     "px": "Location of camera in x-direction",
     "py": "Location of camera in y-direction",
@@ -30,6 +29,7 @@
 parser = argparse.ArgumentParser()
 parser.add_argument("--radiance", action='store_true', help="example settings for computing hemispheric radiance distributions with a sky view camera")
 parser.add_argument("--image", action='store_true', help="example settings for creating visual images with a square camera looking horizontally")
+parser.add_argument("--toa", action='store_true', help="example settings for obtaining top-of-atmosphere radiances")
 parser.add_argument("--sza", type=float, help="solar zenith angle")
 parser.add_argument("--azi", type=float, help="solar azimuth angle")
 parser.add_argument("--name", type=str, default='rte_rrtmgp_input.nc', help="raytracer input file")
@@ -54,17 +54,16 @@
 # add variables: if not available yet
 for var in camera_variables:
     try:
-        cam.createVariable(var,"f8")
+        cam.createVariable(var,"i4" if var in ['cam_type'] else "f8")
     except:
         pass
 
-# example radiance settings
+# example camera radiance settings
 if args['radiance']:
     cam["yaw"][:]   = 0
     cam["pitch"][:] = -90
     cam["roll"][:]  = 0
-    cam["fisheye"][:]= 1
-    cam["f_zoom"][:]= 1
+    cam["cam_type"][:]= 1
     cam["fov"][:] = 80
     cam["px"][:] = 0
     cam["py"][:] = 0
@@ -77,22 +76,27 @@
     cam["yaw"][:] = 0
     cam["pitch"][:] = 0
     cam["roll"][:]  = 0
-    cam["fisheye"][:]= 0
-    cam["f_zoom"][:]= 1
+    cam["cam_type"][:]= 0
     cam["fov"][:] = 80
     cam["px"][:] = 0.
     cam["py"][:] = 0.
     cam["pz"][:] = 500.
     cam["nx"][:] = 256
     cam["ny"][:] = 256
 
+# example toa-radiance settings
+if args['toa']:
+    cam["cam_type"][:]= 2
+    cam["nx"][:] = 256
+    cam["ny"][:] = 256
+
 for var in camera_variables:
     if not args[var] is None:
         cam[var][:] = args[var]
 
 if not args['sza'] is None:
     try:
-        ncf.createVariable('sza','f4',ncf['mu0'].dimensions) 
+        ncf.createVariable('sza','f4',ncf['mu0'].dimensions)
     except:
         pass
     ncf['sza'][:] = np.deg2rad(args['sza'])
@@ -103,8 +107,12 @@
 
 print("Camera settings:")
 for v in camera_variables.keys():
-    print("{:6}{:>8}".format(v, str(cam[v][:])))
-print("{:6}{:>8}".format("sza", str(np.round(np.rad2deg(np.arccos(ncf['mu0'][:].flatten()[0])),1))))
-print("{:6}{:>8}".format("azi", str(np.round(np.rad2deg(ncf['azi'][:].flatten()[0]),1))))
+    if v == 'cam_type':
+        icam = int(cam[v][:])
+        print("{:8}{:>8} ({:s})".format(v, str(icam), ['fisheye','rectangular','TOA radiance'][icam]))
+    else:
+        print("{:8}{:>8}".format(v, str(cam[v][:])))
+print("{:8}{:>8}".format("sza", str(np.round(np.rad2deg(np.arccos(ncf['mu0'][:].flatten()[0])),1))))
+print("{:8}{:>8}".format("azi", str(np.round(np.rad2deg(ncf['azi'][:].flatten()[0]),1))))
 
 ncf.close()

src_kernels_cuda_rt/raytracer_kernels_bw.cu

@@ -196,23 +196,53 @@ namespace
         const Float i = (Float(ij_cam % camera.nx) + rng())/ Float(camera.nx);
         const Float j = (Float(ij_cam / camera.nx) + rng())/ Float(camera.ny);
 
-        if (camera.fisheye)
+
+        if (camera.cam_type == 0)
         {
-            const Float photon_zenith = i * Float(.5) * M_PI / camera.f_zoom;
+            // Fish eye camera
+            const Float photon_zenith = i * Float(.5) * camera.fov / Float(180.) * M_PI;
             const Float photon_azimuth = j * Float(2.) * M_PI;
-            const Vector<Float> dir_tmp = {sin(photon_zenith) * sin(photon_azimuth), sin(photon_zenith) * cos(photon_azimuth), cos(photon_zenith)};
+            const Vector<Float> dir_tmp = {cos(photon_zenith), sin(photon_zenith) * cos(photon_azimuth), sin(photon_zenith) * sin(photon_azimuth)};
 
             photon.direction.x = dot(camera.mx,  dir_tmp);
             photon.direction.y = dot(camera.my,  dir_tmp);
-            photon.direction.z = dot(camera.mz,  dir_tmp) * Float(-1);
+            photon.direction.z = dot(camera.mz,  dir_tmp);
+            photon.position = camera.position + s_min;
         }
-        else
+        else if (camera.cam_type == 1)
         {
             // Rectangular camera based on Villefranque et al. 2019
             photon.direction = normalize(camera.cam_width * (Float(2.)*i-Float(1.0)) + camera.cam_height * (Float(2.)*j-Float(1.0)) + camera.cam_depth);
+            photon.position = camera.position + s_min;
         }
+        else
+        {
+            // Top-of-atmosphere upwelling radiances
+            photon.direction = {Float(0.), Float(0.), Float(-1)};
+            photon.position.x = (Float(ij_cam % camera.nx) + Float(0.5)) * (grid_size.x / camera.nx);
+            photon.position.y = (Float(ij_cam / camera.nx) + Float(0.5)) * (grid_size.y / camera.ny);
+            photon.position.z = z_lev_bg[kbg] - s_min;
+        }
+
+        // if camera starts above domain top, bring photon towards domain top
+        if ((photon.position.z > z_lev_bg[kbg]) and (photon.direction.z < Float(0.)))
+        {
+            Float ds = (photon.position.z - z_lev_bg[kbg]) / photon.direction.z;
 
-        photon.position = camera.position + s_min;
+            photon.position.z = z_lev_bg[kbg] - s_min;
+            photon.position.y += photon.direction.y * ds;
+            photon.position.x += photon.direction.x * ds;
+
+            // Cyclic boundary condition in x.
+            photon.position.x = fmod(photon.position.x, grid_size.x);
+            if (photon.position.x < Float(0.))
+                photon.position.x += grid_size.x;
+
+            // Cyclic boundary condition in y.
+            photon.position.y = fmod(photon.position.y, grid_size.y);
+            if (photon.position.y < Float(0.))
+                photon.position.y += grid_size.y;
+        }
 
         photon.kind = Photon_kind::Direct;
         weight = 1;
@@ -764,23 +794,30 @@ void ray_tracer_kernel_bw(
             Float dist = 0;
             bool reached_cloud = false;
 
-            //const int pix = n * pixels_per_thread + ipix;
             const Float i = (Float(pix % camera.nx) + Float(0.5))/ Float(camera.nx);
             const Float j = (Float(pix / camera.nx) + Float(0.5))/ Float(camera.ny);
 
-            position = camera.position + s_eps;
 
-            if (camera.fisheye)
+            if (camera.cam_type == 0)
             {
                 const Float photon_zenith = i * Float(.5) * M_PI / camera.f_zoom;
                 const Float photon_azimuth = j * Float(2.) * M_PI;
                 const Vector<Float> dir_tmp = {sin(photon_zenith) * sin(photon_azimuth), sin(photon_zenith) * cos(photon_azimuth), cos(photon_zenith)};
 
                 direction = {dot(camera.mx,  dir_tmp), dot(camera.my,  dir_tmp), dot(camera.mz,  dir_tmp) * Float(-1)};
+                position = camera.position + s_eps;
             }
-            else
+            else if (camera.cam_type == 1)
             {
                 direction = normalize(camera.cam_width * (Float(2.)*i-Float(1.0)) + camera.cam_height * (Float(2.)*j-Float(1.0)) + camera.cam_depth);
+                position = camera.position + s_eps;
+            }
+            else
+            {
+                direction = {Float(0.), Float(0.), Float(-1)};
+                position = {i * grid_size.x / camera.nx,
+                            j * grid_size.y / camera.ny,
+                            grid_size.z - 2*s_eps};
             }
 
             // first bring photon to top of dynamical domain

src_test/Radiation_solver_bw.cu

@@ -630,6 +630,8 @@ void Radiation_solver_longwave::solve_gpu(
     */
 }
 
+// color matchin functions from https://jcgt.org/published/0002/02/01/paper.pdf
+
 
 Float get_x(const Float wv)
 {