DATA_SECTION
  init_int nbins
  init_vector length_bins(1,nbins);
  int ilen
  int il2
  number dev

PARAMETER_SECTION
  init_bounded_number af1(0.5656024,0.56560241,1)   // Female ln(a) for ln(a) *len^b  growth-increment
  init_bounded_number bf1(0.913266,0.9132661,-8)    // Female ln(a) for ln(a) *len^b  growth-increment
  init_bounded_number beta(0.75000,0.75001,-2)      // Growth beta
  number alpha;
  matrix len_len(1,nbins,1,nbins);
  vector mean_length(1,nbins);
  objective_function_value objf

PROCEDURE_SECTION
  for(ilen=1;ilen<=nbins;ilen++)
  {
    mean_length(ilen)= mfexp(af1)* pow(length_bins(ilen),bf1);
    alpha = (mean_length(ilen)-(length_bins(ilen)-2.5))/beta;
    for(il2=ilen;il2<=nbins;il2++)
    {
      dev = length_bins(il2)+2.5-length_bins(ilen);

      // This line works
      // len_len(ilen,il2) = pow(dev,(alpha-1.)) *   exp(-dev/beta);

      // This line fails to produce the appropriate matrix
      len_len(ilen,il2) = pow(dev,(alpha-1.)) * mfexp(-dev/beta);

      // This line also works
      // len_len(ilen,il2) = mfexp(log(dev)*(alpha-1.) - dev/beta);
    }  
    len_len(ilen) /= sum(len_len(ilen));
  }
  cout <<len_len<<endl;exit(1);
