#include "stdio.h"
#include "math.h"
#include "string.h"

//  This program calculates the Analog Devices 9830A DDS frequency divisor
//  values for programming the FEI FE-5650A and FE-5680A rubidium oscillator
//  units.
//
//  This program was written by Mark S Sims  May 2008.	It is relased to the
//  public domain.  Use it in peace.  Plase share with others what you do
//  with it.
//

// These are the reference freq and divisors returned by the FEI "S" command
// They vary from unit to unit.  Change these to whatever your unit returns (or
// modify this program to interrogate the unit with the "S" command.  While you
// are at it,  modify this program to send the divisor word to the oscillator.
//
#define R  50255055.011982
#define F  0x2ABB5050L

#define FF ((double) 8388608.0)      // the freq the module was calibrated for
#define M  ((double) 4294967296.0)   // 2^32

double r=R;    // working reference frequency
double f=FF;   // working DDS divisor word


void lll(double x)
{  // print divisor as a 64 bit hex value
char s[20];
int i;

   for(i=0; i<16; i++) s[i] = 0;
   i = 0;

   x = x * M * M;
   while(x>0.99) {
      s[i++] = fmod(x, 16);
      if(i > 16) break;
      x = x / 16.0;
   }

   for(i=15; i>=0; i--) {
      printf("%X", s[i]);
      if(i == 8) printf(":");
   }
   printf("\n\n");
}

main(int argc, char *argv[])
{
double d;
double x;

   if(argc > 1) {  // get desired freq from command line
      sscanf(argv[1], "%lf", &f);
      if(argc > 2) {  // if more than one command line arg,  freq is in MHz or KHz
	 if(toupper(argv[2][0]) == 'K') f *= (double) 1000.0;
	 else if(toupper(argv[2][0]) == 'H') f *= (double) 1.0;
	 else f *= (double) 1000000.0;
      }
      else if(strchr(argv[1], 'K')) f *= (double) 1000.0;
      else if(strchr(argv[1], 'k')) f *= (double) 1000.0;
      else if(strchr(argv[1], 'M')) f *= (double) 1000000.0;
      else if(strchr(argv[1], 'm')) f *= (double) 1000000.0;
   }

   // print input data values
   printf("\nR=%19.10lf  F=%08lX   (f=%lf)\n\n",  R,F,FF);

   // calculate what freq the R= and F= values returned from the
   // oscillator would produce
   d = ((double) F) / M;
   d *= r;
   printf("R*F/(2^32)=%19.10lf  ref_scale=%.14lf\n", d, FF/d);

   // scale the R= reference freq by that value to get the actual
   // oscillator reference freq
   r = r * FF / d;
   printf("actual ref=%.10lf\n\n", r);

   // show what the 64 bit DDS divisor should be
   printf("freq=%.10lf\n\n", f);
   d = f / r;
   printf("freq/ref=%.19lf\n", d);

   printf("full divisor=");
   lll(d);


   // calculate the closest 32 bit divisor that the AD9830A DDS chip
   // will actually use
   d *= M;
   d = (double) (unsigned long) (d+0.5);


   // now show the actual generated frequency and the two closest ones
   x = r * (long) (d-1);
   x /= M;
   printf("freq(%08lX) = %.10lf  err=%13.10lf (%.3lg)\n", (long)d-1, x,  x-f, (x-f)/f);

   x = r * (long) d;
   x /= M;
   printf("freq(%08lX) = %.10lf  err=%13.10lf (%.3lg)\n", (long)d, x,  x-f, (x-f)/f);

   x = r * (long) (d+1);
   x /= M;
   printf("freq(%08lX) = %.10lf  err=%13.10lf (%.3lg)\n", (long)d+1, x,  x-f, (x-f)/f);

   printf("\n");
   printf("nearest divisor=%08lX\n", (unsigned long) d);
}