/* Last edited on 2024-10-21 05:36:10 by stolfi */

/* Version of {multifok_rayset_choose_tilts_and_weights} that
  is simply a digital disk of radius {HS}. */

void multifok_rayset_choose_tilts_and_weights
  ( int32_t HR, 
    int32_t *NR_P,
    r2_t **tRay_P, 
    double **wRay_P,
    bool_t verbose
  )
  {
    /* This implementation chooses ray tilts in a regular orthogonal grid pattern with round border. */
    
    demand(HR >= 0, "invalid {HR}");

    int32_t NR_max = (2*HR+1)*(2*HR+1);
    
    /* Allocate {tRay,wRay} for full square grid: */
    r2_t *tRay = talloc(NR_max, r2_t);
    double *wRay = talloc(NR_max, double);
    
    double sigma = 0.50; /* Deviation of weight function. */
    double r2max = 1.05*(1.0 + 1.0/((double)HR)); /* Trimming limit. */
    
    int32_t NR = 0; /* Rays generated so far. */
    
    /* Central ray is vertical: */
    tRay[NR] = (r2_t){{ 0, 0 }};
    wRay[NR] = 1.0;
    NR++;
    
    if (HR > 0)
      { /* Fill rest of {tRay,wRay} with digital disk only, minus vertical ray: */
        for (int32_t iy = -HR; iy <= +HR; iy++)
          for (int32_t ix = -HR; ix <= +HR; ix++)
            { if ((ix != 0) || (iy != 0))
                { double x = ((double)ix)/HR;
                  double y = ((double)iy)/HR;
                  double r2 = x*x + y*y;
                  if (r2 <= r2max)
                    { assert(NR < NR_max);
                      tRay[NR] = (r2_t){{ x, y }};
                      wRay[NR] = exp(-(x*x + y*y)/(sigma*sigma)/2);
                      NR++;
                    }
                }
            }
        assert(NR >= 5);

        /* Normalize the weights to unit sum, and the titls to unit RMS radius: */
        double sum_w = 0;
        double sum_w_r2 = 0;
        for (int32_t ir = 0; ir < NR; ir++)
          { sum_w += wRay[ir];
            double r2 = r2_norm_sqr(&(tRay[ir])); 
            sum_w_r2 += wRay[ir]*r2;
           }
        assert(sum_w >= 1.0);
        double r2_avg = sqrt(sum_w_r2/sum_w);
        assert(r2_avg > 0.0);
        for (int32_t ir = 0; ir < NR; ir++) 
          { wRay[ir] /= sum_w; 
            tRay[ir].c[0] /= r2_avg;
            tRay[ir].c[1] /= r2_avg;
          }
      }
    
    /* Resize the arrays: */
    tRay = realloc(tRay, NR*sizeof(r2_t));
    wRay = realloc(wRay, NR*sizeof(double));

    if (verbose)
      { fprintf(stderr, "ray tilts and weights:\n");
        for (int32_t ir = 0; ir < NR; ir++)
          { fprintf(stderr, "      ray %4d tRay = ", ir);
            r2_gen_print(stderr, &(tRay[ir]), "%+9.6f", "( ", " ", " )");
            fprintf(stderr, " wRay = %12.10f\n", wRay[ir]);
          }
      }

    /* Return results: */
    (*NR_P) = NR;
    (*tRay_P) = tRay;
    (*wRay_P) = wRay;
  }