#include <cmacs.h>
#include <image.h>
#include <imagemac.h>

IMAGE *pic,*picr;   int sh;   

main(argc,argv)
int argc; char *argv[];
{  int     i,j,prl,prh,pcl,pch,frl,fcl,pb,rv,rh,rt,ii,jj,rvb,rhb,pixelmax;
   float   lat,lng,lnglo,lnghi,latlo,lathi,x,y,z,r,dlat,dlng;
   int line[1024];  char *lineo;

   rv=40;
   rh=rv*0.75*512/480;
   rvb = rv/5;  rhb = rh/5;

   pic = imgopen("EARTH.img",0);
   pixelmax = (1 << pixbits(pic)) - 1;

     printf("%d rows x %d cols x %d bits\n",rows(pic),cols(pic),pixbits(pic));
     prl = rstart(pic);  prh = rend(pic);
     pcl = cstart(pic);  pch = cend(pic);
     pb = pixbits(pic);
     picr = imgcreat("synch.img",0655,pb,0,479,0,511);
     for (jj=rh; jj <= 511-rhb; ++jj)
      { ii = (479.0-rv) + (rvb - 479.0 + rv)*(jj-rh)/(511.0-rhb-rh);
        setpixel(picr,ii,jj,pixelmax); }
     for (ii = -rv; ii <= rv; ++ii)
       {if ((ABS(ii) % 10) == 0) {printf("ii %d ",ii); fflush(stdout);}
       y = ((double) ii)/((double) rv);
       rt = SQRT(1.0-y*y)*rh;
       for (jj = -rt; jj <= rt; ++jj)
          {x = ((double) jj)/((double) rh);   r = SQRT(x*x+y*y);
           z = SQRT(1.0-r*r);
           lng = ATAN2(x,y);   lat = ATAN2(z,r);
           lng = lng + 1.4*PI/2;
           if (lng>PI) lng -= TWOPI;
           j = (lng+PI)*(pch-pcl)/TWOPI + pcl;
           i = (-lat+PI/2)*(prh-prl)/PI + prl;
           j = MIN(pch,MAX(pcl,j));
           i = MIN(prh,MAX(prl,i));
           setpixel(picr,ii+479-rv,jj+rh,pixel(pic,i,j));
           }
         }

     for (ii = -rvb; ii <= rvb; ++ii)
       {y = ((double) ii)/((double) rvb);
       rt = SQRT(1.0-y*y)*rhb;
       for (jj = -rt; jj <= rt; ++jj)
          {setpixel(picr,ii+rvb,jj+511-rhb,pixelmax);}
        }
       imgclose(picr);  imgclose(pic);

}