// RAMAN-NATH EQUATIONs
// i/pi*dS_j/dz = j*(th-j)*S_j+a*(S_(j+1)+S_(j-1))
// sinusoidal grating diffraction:  power in j-th mode (input in 0-th)
// rn(j=order, s=z/Ltalbot, a=Ltalbot/Lpendell., th=angle(1=bragg))
// (Ltalb=period^2*n/wvlgth0,  Lpend=2*wlgth*sqrtepsilon/deltaepsilon,
//  braggangle=wvlgth0/(2*period) )                       (j may be array)
// input: row vectors, output: column vectors or matrix(j,s)
                                       // (accuracy of the order of 0.001)
clear
function yy=rana(j,s,a,th);
  if length(a)>1, 
    for h=1:length(a), yy(:,h)=rn(j,s,a(h),th); end
  else if length(th)>1, 
    for h=1:length(th), yy(:,h)=rn(j,s,a,th(h)); end
    else yy=rn(j,s,a,th);
    end, end
    endfunction
     
     function xx = rn(j,s,a,th);
     if or((j-round(j)) ~= 0), 
        error('order must be integer'), 
     end,
     les = length(s);
     if a==0, xx = zeros(length(j),les);                           // no grating
        if min(j)<1, 
           xx(-min(j)+1,:) = ones(1,les); 
        end,
     else 
        if max(s)<0.033, xx = besselj(j',2*%pi*s*a).^2;     // Raman-Nath limit
          xx=xx';
        else       
           n0=2;                            // increase n0 for higher precision
           f = s*a;                                                  // f=z/Lp
              if max(f)*s<0.1, n=fix(max(f)*10)+n0;  
              else n = fix(2*sqrt(a))+n0;              // 2*n+1 = matrix order   
              end,                             // (< 0.001 power in orders >n)
           Y = ones(2*n+1,1);
           C = zeros(2*n+1,1);
           C(n+1) = 1;                  // initial condition: all zero-th mode
           arr = %pi*(-n:n).*(th-(-n:n))+2*%pi*a;
           sub = ones(1,2*n)*%pi*a;
           H = diag(arr) + diag(sub,1) + diag(sub,-1);
           [U,D] = spec(H);                                       // eigenmodes
           for k=1:les,
            S(:,k) = U * (exp(-%i*s(k)*D * Y) .* (inv(U) * C));   // propagation
              //if (j<-n)or(j>n), xx = zeros(length(j),les);    // if orders all >n
              //else 
              jj = max(min(j),-n):min(max(j),n);
                 SS = (abs(S(n+1+jj,k))).^2;
                 z1 = zeros(-min(j)-n,1);         z2 = zeros(max(j)-n,1);
                 xx(:,k) = [z1; SS; z2];
              //end,
            end,
        end,
       if length(j)==1, xx = xx';  end
   end
   endfunction
