#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#define PI (4*atan(1.0))
#define DEBUG 0

void findrates(double *azrate, double *elrate, double latrad, double elrad, double azrad);

void rateusage(char *a, double latitude, double minelev, double azmaxrate, double elmaxrate) {
  printf("Usage: %s [-r <rateDegPerSec>] -a <azimDeg> -e <elevDeg> -z <zaDeg>\n",a);
  printf("   -M -l <latDeg:default=%.4f> -m <minSmoothVel(deg/min)> -grid -s <stepDeg> -VLA\n", latitude);
  printf("   -E <minElevDeg:default=%.2f>  --nofile  -findaz\n", minelev);
  printf("   --azmaxrate <default=%.0f deg/min>  --elmaxrate <default=%.0f deg/min>\n",azmaxrate,elmaxrate);
  printf("Usage1: -r -a   to compute maximum elevation given rate r\n");
  printf("Usage2: -a -e   to compute az and el tracking rates\n");
  printf("Usage3: -grid   to compute az and el tracking rates on a grid\n");
  printf("Usage4: -R      to compute az and el tracking rates from a random distribution of az & el\n");
  printf("Usage5: -findaz to compute az,el combinations with zero az rate\n");
  printf("  -M selects the latitude of Mauna Kea\n");
  printf("  -VLA selects the latitude of the VLA\n");
  exit(0);
}

int main(int argc, char *argv[]) {
  char outfile[100];
  double rate, azrate, elrate, msv, azminrate, elminrate, azmaxrate, dec;
  double az,z,latitude,el,zrad, azratedegmin, elratedegmin, minelev ;
  double azrad, elrad, latrad, azratedeg, elratedeg, angstep, elmaxrate;
  double prevaz, prevazrate, prevazrad, sky, sqdeg;
  int i, findrate = 1, azpoints, azonly=0, elonly=0, found;
  int debug = DEBUG;
  FILE *f;
  int dogrid = 0, success, points, total, toofast;
  int nofile = 0;
  int findaz = 0;
  az = 0;
  latitude = 38.43294; 
  el = latitude;
  msv =  0.06; /* deg/min */
  azmaxrate = 36; /* deg/min */
  elmaxrate = 18; /* deg/min */
  angstep = 0.1;
  minelev = 5; /* GBT */

  for (i=1; i<argc; i++) {
    if (strstr(argv[i],"-d")!=NULL) {
      debug = 1;
      continue;
    }
    if (strstr(argv[i],"-nofile")!=NULL) {
      nofile = 1;
      continue;
    }
    if (strstr(argv[i],"-findaz")!=NULL) {
      findaz = 1;
      continue;
    }
    if (strstr(argv[i],"-grid")!=NULL) {
      dogrid = 1;
      continue;
    }
    if (strstr(argv[i],"-VLA")!=NULL) {
      latitude = 34.02777; /* deg = VLA */
      continue;
    }
    if (strstr(argv[i],"-R")!=NULL) {
      dogrid = 2;
      continue;
    }
    if (strstr(argv[i],"-r")!=NULL) {
      if (++i < argc) {
	sscanf(argv[i],"%lf",&rate);
	findrate = 0;
	continue;
      } else {
	printf("-r requires an argument\n");
	exit(0);
      }
      continue;
    }
    if (strstr(argv[i],"-h")!=NULL) {
      rateusage(argv[0],latitude,minelev,azmaxrate,elmaxrate);
    }
    if (strstr(argv[i],"--azonly")!=NULL) {
      azonly = 1;
      continue;
    }
    if (strstr(argv[i],"--elonly")!=NULL) {

      continue;
    }
    if (strstr(argv[i],"--azmaxrate")!=NULL) {
      if (++i < argc) {
	sscanf(argv[i],"%lf",&azmaxrate);
	continue;
      } else {
	printf("--azmaxrate requires an argument\n");
	exit(0);
      }
      continue;
    }
    if (strstr(argv[i],"--elmaxrate")!=NULL) {
      if (++i < argc) {
	sscanf(argv[i],"%lf",&elmaxrate);
	continue;
      } else {
	printf("--elmaxrate requires an argument\n");
	exit(0);
      }
      continue;
    }
    if (strstr(argv[i],"-E")!=NULL) {
      if (++i < argc) {
	sscanf(argv[i],"%lf",&minelev);
	continue;
      } else {
	printf("-E requires an argument\n");
	exit(0);
      }
      continue;
    }
    if (strstr(argv[i],"-a")!=NULL) {
      if (++i < argc) {
	sscanf(argv[i],"%lf",&az);
	continue;
      } else {
	printf("-a requires an argument\n");
	exit(0);
      }
    }
    if (strstr(argv[i],"-e")!=NULL) {
      if (++i < argc) {
	sscanf(argv[i],"%lf",&el);
	z = 90-el;
	continue;
      } else {
	printf("-e requires an argument\n");
	exit(0);
      }
    }
    if (strstr(argv[i],"-z")!=NULL) {
      if (++i < argc) {
	sscanf(argv[i],"%lf",&z);
	continue;
      } else {
	printf("-z requires an argument\n");
	exit(0);
      }
    }
    if (strstr(argv[i],"-M")!=NULL) {
      latitude = +19.82547; /* old CSO value */
      continue;
    }
    if (strstr(argv[i],"-l")!=NULL) {
      if (++i < argc) {
	sscanf(argv[i],"%lf",&latitude);
	continue;
      } else {
	printf("-l requires an argument\n");
	exit(0);
      }
    }
    if (strstr(argv[i],"-m")!=NULL) {
      if (++i < argc) {
	sscanf(argv[i],"%lf",&msv);
	continue;
      } else {
	printf("-m requires an argument\n");
	exit(0);
      }
    }
    if (strstr(argv[i],"-s")!=NULL) {
      if (++i < argc) {
	sscanf(argv[i],"%lf",&angstep);
	continue;
      } else {
	printf("-s requires an argument\n");
	exit(0);
      }
    }
    printf("Unrecognized argument = %s\n",argv[i]);
    exit(0);
  }
  sky = 18828.76;  /* 2*pi*(1-sin(5*pi/180)) * pow(180/pi,2) */
  latrad = latitude * PI/180.0;
  if (dogrid>0) {
    if (nofile==0) {
      strcpy(outfile,"azelrate.out");
      f = fopen(outfile,"w");
      fprintf(f,"# Az El AzRate(deg/min) ElRate(deg/min)  latitude %f\n",latitude);
    }
    azminrate = elminrate = msv;
    success = 0;
    points = 0;
    total = toofast = 0;
    el = minelev;
    while (el <= 90) {
      az = 0;
      if (debug) {
	printf("El = %.1f ",el);
      }
      azpoints = 0;
      while (az < 360) {
	elrad = el*PI/180.0;
	az += angstep/cos(elrad);
	azrad = az*PI/180.0;
	if (dogrid==2) {
	  elrad = (5 + 85.0*rand()/2147483647)*PI/180.;
	  azrad = 2*PI*(1.0*rand()/2147483647);
	}
	findrates(&azrate,&elrate, latrad, elrad, azrad);
	azratedegmin = azrate*180*60/PI;
	elratedegmin = elrate*180*60/PI;
	if (nofile==0) {
	  fprintf(f,"%.1f %.1f %f %f\n",az,el,azratedegmin,elratedegmin);
	}
	if (fabs(azratedegmin) < azmaxrate && fabs(elratedegmin) < elmaxrate) {
	  points++;
	} else {
	  toofast++;
	  /*	  printf("Too fast = %d = %f\n",toofast,azratedegmin);*/
	}
	total++;
	if (azonly == 1) {
	  if (fabs(azratedegmin)>=azminrate) {
	    success++;
	  }
	} else if (elonly == 1) {
	  if (fabs(elratedegmin)>=elminrate) {
	    success++;
	  }
	} else {
	  if (fabs(azratedegmin)>=azminrate && fabs(elratedegmin) >= elminrate) {
	    success++;
	  }
	}
	azpoints++;
      }
      if (debug) {
	printf("(%3d) ",azpoints);
      }
      el += angstep;
    }
    if (debug) {
      printf("\n");
    }
    printf("For MSV = %.3f deg/min, %d/%d = %.3f%% success\n",msv,success,total,success*100.0/total);
    sqdeg = sky*toofast/total;
    /* 2*PI*h*pow(180/PI,2) = sqdeg
       r = sin(el);  h=1-r;  h=1-sin(el)
       2*PI*(1-sin(el))*pow(180/PI,2) = sqdeg
       1-sin(el) = sqdeg/(2*PI*pow(180/PI,2))
       sin(el) = 1-(sqdeg/(2*PI*pow(180/PI,2)))
       el = asin(1-(sqdeg/(2*PI*pow(180/PI,2)))
    */
    el = 180*asin(1-(sqdeg/(2*PI*pow(180/PI,2))))/PI;
    fprintf(stderr,"  toofast=%d=%.5f%%, times %.2f sq deg w/el>5 = %f sq deg\n",toofast,
	    (toofast*100.0)/total,sky,   sqdeg);
    fprintf(stderr," = (za>%.4f) = el<%.4f\n",sqrt(sqdeg/PI),el);
    if (nofile==0) {
      fclose(f);
      printf("Results written to: %s\n",outfile);
    }
    exit(0);
  }
  if (findaz > 0) {
    if (nofile==0) {
      f = fopen("azelrate.zero","w");
    }
    az = 73.0;
    do {
      azrad = az*PI/180.;
      found = 0;
      printf("Checking az=%f\n",az);
      for (el=minelev; el<90; el+=angstep) {
	elrad = el*PI/180.;
	prevazrate = azrate;
	prevazrad = azrad;
	prevaz = az;
	findrates(&azrate,&elrate, latrad, elrad, azrad);
	if (el>minelev) {
	  if ((azrate <= 0 && prevazrate > 0) || (azrate >= 0 && prevazrate < 0)) {
	    /* compute dec */
            dec = (180./PI)*asin(sin(latrad)*sin(elrad) + cos(latrad)*cos(elrad)*cos(0.5*(prevazrad+azrad)));
	    printf("%f %f %+f\n",0.5*(prevaz+az),el,dec);
	    if (nofile==0) {
	      fprintf(f,"%f %f %+f\n",0.5*(prevaz+az),el,dec);
	    }
	    found = 1;
	    break;
	  }
	}
      }
      az += angstep;
    } while (found);
    az = 73.0-angstep;
    do {
      azrad = az*PI/180.;
      found = 0;
      printf("Checking az=%f\n",az);
      for (el=minelev; el<90; el+=angstep) {
	elrad = el*PI/180.;
	prevazrate = azrate;
	prevaz = az;
	prevazrad = azrad;
	findrates(&azrate,&elrate, latrad, elrad, azrad);
	if (el>minelev) {
	  if ((azrate <= 0 && prevazrate > 0) || (azrate >= 0 && prevazrate < 0)) {
            dec = (180./PI)*asin(sin(latrad)*sin(elrad) + cos(latrad)*cos(elrad)*cos(0.5*(prevazrad+azrad)));
	    printf("%f %f %+f\n",0.5*(prevaz+az),el,dec);
	    if (nofile==0) {
	      fprintf(f,"%f %f %+f\n",0.5*(prevaz+az),el,dec);
	    }
	    found = 1;
	    break;
	  }
	}
      }
      az -= angstep;
    } while (found);
    if (dogrid == 2) {
      fprintf(stderr,"Warning: This calculation used random distribution of az,el, which\n");
      fprintf(stderr,"         over-emphasizes the area near the zenith.  Use -grid instead.\n");
    }
    fclose(f);
    exit(0);
  }

  /* if we get this far, then we are just asking for a single position */
  zrad = z*PI/180.0;
  azrad = az*PI/180.0;
  elrad = el*PI/180.0;
  if (findrate) {
    findrates(&azrate,&elrate, latrad, elrad, azrad);
    azratedeg = azrate*180/PI; /* deg/sec */
    elratedeg = elrate*180/PI;
    azratedegmin = azratedeg*60;
    elratedegmin = elratedeg*60;
    printf("For az=%f, el=%f\n", az, el);
    printf("azrate = %.4f deg/sec = %f deg/min\n",azratedeg, azratedegmin);
    printf("elrate = %.4f deg/sec = %f deg/min\n",elratedeg, elratedegmin);
  } else {
    rate = 1.32; /* deg/sec */
    rate *= 60;  /* deg/min */
    rate *= 4;
    rate += sin(latitude);
    rate *= cos(latitude)*cos(az); /* now equal to cot(z) */
    z = atan(1.0/rate);
    z *= 180/PI;
    printf("ZA = %f deg\n",z);
  }
  return(0);
}

void findrates(double *azrate, double *elrate, double latrad, double elrad, double azrad) {
    /* compute azimuth and elevation rate  in radians per second at the specified az and el. */
    /* use eq from http://www.astro.caltech.edu/~mcs/CBI/pointing/ */
    /* replace "cos(decrad)*cos(prad)" with 
     *         "sin(latrad)*cos(elrad)-cos(latrad)*sin(elrad)*cos(azrad)"
     */
  double sidereal = (2*PI/86400)*1.002737811906;
#if 1
  *azrate =  sidereal*((sin(latrad)*cos(elrad)-cos(latrad)*sin(elrad)*cos(azrad))/cos(elrad));
#else
  *azrate = sin(latrad) - sin(
#endif
  *elrate =  sidereal*(cos(latrad)*sin(azrad));
}