// fra0.c  fraunhofer rings and low elevation SSPS rectenna.
// cc -o fra0 fra0.c -lm
//
// This computes the path from an -5° SSPS to a rectenna at 46° N 
// and 30° west.  This gives us the antenna elevation.  We then 
// assume a circular aperture, and a rectenna expanded by 1/sin(elevation)
// to catch diagonal beams.  
//
// Skolnik says 0.7dB loss at 2.5 GHz and 25 degrees elevation, 
//   thus an 85% efficiency multiplier.  If the rectenna is 85% 
//   efficient, and we have 91% beam capture efficiency and 95%
//   electrical grid efficiency, the broadcast to load efficiency
//   is 69% .  That turns 7 GW broadcast into 4.8 GW to the customer.

#define	 STEP	0.001		// m q step
#define  DMAX	2113.1000	// m   rectenna diameter
#define	 APER	5000.0000	// m   ssps aperture diameter

#define  FREQ	2.450E+09	// hz  frequency
#define	 C	2.998E+08	// m/s lightspeed
#define	 GEO	42164000	// m GEO radius
#define  RE	6371000  	// m Earth radius
#define  SOUTH	5.00    	// degrees SSPS offset latitude south 
#define  LAT	46.1920  	// degrees rectenna latitude
#define	 LON	30.0000		// degrees rectenna longitude west

#include <math.h>
#include <stdio.h>

int main() {
   double pi   = 4.0*atan(1.0)			 ;//       pi
   double d2r  = atan(1.0)/45.0			 ;//       degrees to radians
   double r2d  = 45.0/atan(1.0)			 ;//       radians to degrees
   double dNS  = SOUTH+LAT			 ;// deg   northsouth
   double rNS  = d2r*dNS	         	 ;// rad   northsouth
   double wl   = (C/FREQ)			 ;// m     wavelength
   double rEW  = d2r*LON			 ;// rad   eastwest
   double rtot = acos(cos(rNS)*cos(rEW))	 ;// rad   total angle
   double dtot = r2d*rtot			 ;// deg   total angle
   double zz   = RE*sin(rtot)			 ;// m     Z
   double xg   = GEO-RE*cos(rtot)		 ;// m     X to geo
   double dist = sqrt( zz*zz + xg*xg )		 ;// m     distance 
   double rel  = atan(xg/zz)-rtot		 ;// rad   elevation
   double del  = r2d*rel			 ;// deg   elevation
   double d    = APER				 ;// m     aperture diameter
   double xsca = pi * d / (wl*dist)		 ;// 1/m   scaling argument

   double qs   = STEP				 ;// m     integration step
   double q    					 ;// m     center of annulus

   double i0p  = 2.0*qs/xsca*xsca		 ;// 1/m2  power scaling
   double pow  = 0.0				 ;//       total power
   double bs                                     ;//       bessel function
   double btest                                  ;//
   double rmax = 0.5*DMAX			 ;// m     rectenna diameter

   for( q = 0.5*qs ; q < DMAX ; q += STEP ) {
      bs   = j1(xsca*q)       			 ;// m2	   additional area
      pow += bs * bs / q			 ;// 
   }
   pow *= i0p 					 ;//
   double lost = 1.0 - pow                       ;//
   double sin1 = 1/sin(rel)                      ;//
   double sin2 = DMAX/sin(rel)                   ;//

   printf( "%10.3e    bessel\n"                  , bs    ) ;
   printf( "%10.7f    fractional power\n"        , pow   ) ;
   printf( "%10.7f    lost power\n"              , lost  ) ;
   printf( "%9.1f m   Rectenna diameter\n"       , DMAX  ) ;
   printf( "%9.1f m   SSPS diameter\n"           , APER  ) ;
   printf( "%9.4f°    SSPS south\n"              , SOUTH ) ;
   printf( "%9.4f°    Rectenna latitude\n"       , LAT   ) ;
   printf( "%9.4f°    Rectenna longitude\n"      , LON   ) ;
   printf( "%9.4f°    Rectenna total angle\n"    , dtot  ) ;
   printf( "%9.4f°    Rectenna elevation\n"      , del   ) ;
   printf( "%9.5f     Rectenna 1/sin\n"          , sin1  ) ;
   printf( "%9.1f     Rectenna D/sin\n"          , sin2  ) ;
   printf( "%10.3e m  distance\n"                , dist  ) ;
   return(0);
}