DATA _NULL_;
   length foot $50; 
    CALL SYMPUT ('dayname', SYSGET("EMSDAY"));
    CALL SYMPUT ('casename', SYSGET("CASE"));
/* PERSONALIZE: Change these variables if you would like to change the location of the Legend  */ 
    CALL SYMPUT ('minlat',"25");
    CALL SYMPUT ('minlong',"70");
    CALL SYMPUT ('polid', COMPRESS(UPCASE(SCAN("&sysparm",1,";"))));

/* PERSONALIZE: These variables alter the scale of the circles drawn on the map. Increasing the size of
      values up will increase the size of the circle for a given mass of emissions. */ 
     if COMPRESS(UPCASE(SCAN("&sysparm",1,";"))) = "NOX" then scale = 14000; 
     if COMPRESS(UPCASE(SCAN("&sysparm",1,";"))) = "SO2" then scale = 20000; 
     if COMPRESS(UPCASE(SCAN("&sysparm",1,";"))) = "ROG" then scale = 2000; 
     if scale = . then scale = 14000; 
    CALL SYMPUT ('scale',put(scale,10.));

/* PERSONALIZE: Change these variables to get different output PS=Postscript GIF=GIF image  */ 
    
    device = COMPRESS(UPCASE(SCAN("&sysparm",2,";")))  ; 
    PUT "Device1=" device; 
    if device ne "PS" and device ne "GIF" then device = "GIF" ;    
    CALL SYMPUT ('device',PUT(device,$7.)); 
    PUT "Device2=" device; 

    now = today(); 
    foot = "Created with EMS-2001 on " || trim(compress(PUT(now,monname9.))) || " "  || 
                                          put(now,DAY2.) || "," || Put(now,year4.); 
    PUT "FOOT IS:" foot; 
    CALL SYMPUT ('foot',put(foot,$50.));
    STOP;
  RUN;


%MACRO defgop(dev); 
/* PERSONALIZE: Change these variables to get different output PS=Postscript GIF=GIF image  */ 

   GOPTIONS RESET = goptions;

   GOPTIONS GSFMODE = append
             GSFLEN  = 80
             ROTATE = LANDSCAPE
             ftext=swiss
             GACCESS = gsasfile
             cback=white
             ctext=black
             colors= (black)
             htitle = 4
             htext = 1;


   %IF "&dev" = "PS" %THEN %DO;
      GOPTIONS GPROLOG = '25210d0a'x
               DEVICE = PSCOLOR; 
               FILENAME gsasfile "circlemap.&polid..&dayname..&casename..ps" ;
   %END;
   
   %ELSE %IF "&dev" = "GIF" %THEN %DO;
      GOPTIONS DEVICE = IMGGIF; 
               FILENAME gsasfile "circlemap.&polid..&dayname..&casename..gif" ;
   %END;

   %ELSE %DO;
      GOPTIONS DEVICE = IMGGIF; 
               FILENAME gsasfile "circlemap.&polid..&dayname..&casename..gif" ;

   %END;

%MEND defgop; 

%defgop(&device); 

DATA xy_cord;
     INFILE '$EMS_GRD/project_xy.pt' DLM=',';
     LENGTH fcid $15. stkid $12.;
     INPUT stid cyid fcid $ stkid $ xcoord ycoord lat long;
    /* convert degrees to radians */
    long = -long; 
    x=atan(1)/45 * long; 
    y=atan(1)/45 * lat;
run;

proc sort data=xy_cord; 
 by stid cyid fcid; 
run; 

proc summary data=xy_cord; 
by stid cyid fcid; 
var x y lat long;
output out=coord median=; 
run; 

data ptemis(keep=stid cyid fcid aceeton); 
 set ems_run.ptemis; 
aceeton = aceekg / 907.185; 
if polid ne "&polid" then delete; 
run; 

proc sort data=ptemis; 
by stid cyid fcid; 
run; 

proc summary data=ptemis; 
by stid cyid fcid; 
var aceeton; 
output out=ptemis sum=; 
run; 

data facil; 
 merge coord ptemis; 
by stid cyid fcid; 
if aceeton le 5 or x = . or y = . then delete; 
run; 

data maxs; 
long = &minlong;
 do i = 1 to 4; 

   /*PERSONALIZE: This command calculates the size of the circles in the legend*/ 
   aceeton = int((&scale*i/120)/10) * 10 ; 
   lat = &minlat +2 + (10-(2*i)); 
    x=atan(1)/45 * long;
    y=atan(1)/45 * lat;
  output maxs; 
 end; 
STOP; 
run; 

 data annoleg;
    set maxs;
    retain xsys ysys '2' flag 1 when 'a';
    length text $25 color function $8;

    /* Create an observation to label the service area */
    function='label';
    style='swiss';
    position='2';
    color='black';
    text=trim(left(aceeton))||' Tons/Day';
    size=.75;
    output;

    d2r=3.1415926/180;
    r=&scale;

    /* point for the circle to be drawn around */
    xcen=long;
    ycen=lat;
    radius = aceeton  ** .5;
    /* calculate the points to draw a circle */
    do degree=0 to 360 by 5;
       /* begin a new circle */
       if degree=0 then do;
             /* the circles will be solid yellow */
         function='move';
             line=1;
             size=.45;
             color='black';
           end;
           /* continue drawing the circle */
           else do;
         function='draw';
             /* outline the circle with black */
             color='black';
           end;

           /* calculate a point along the circle */
       y=arsin(cos(degree*d2r)*sin(aceeton/R)*cos(ycen*d2r)+
         cos(aceeton/R)*sin(ycen*d2r))/d2r;
       x=xcen+arsin(sin(degree*d2r)*sin(aceeton/R)/cos(y*d2r))/d2r;

           /* convert degrees to radians */
       x=atan(1)/45*x;
       y=atan(1)/45*y;
           output;
    end;
  run;

data anno;
    set facil;
    retain xsys ysys '2' flag 1 when 'a';
    length text $25 color function $8;

    /* Create an observation to label the service area */
    d2r=3.1415926/180;
    /* radius of the earth in milesi, if you want radius in miles */
    /*  r=3958.739565;*/
    r=&scale; 

    /* point for the circle to be drawn around */
    xcen=long;
    ycen=lat;
    radius = aceeton  ** .5; 
    /* calculate the points to draw a circle */
    do degree=0 to 360 by 5;  
       /* begin a new circle */
       if degree=0 then do;
             /* the circles will be solid yellow */
         function='move';
             line=1;   
             size=.45;
             color='black';
           end;
           /* continue drawing the circle */
           else do;
         function='draw';
             /* outline the circle with black */
             color='black';
           end;

           /* calculate a point along the circle */
       y=arsin(cos(degree*d2r)*sin(aceeton/R)*cos(ycen*d2r)+
         cos(aceeton/R)*sin(ycen*d2r))/d2r;
       x=xcen+arsin(sin(degree*d2r)*sin(aceeton/R)/cos(y*d2r))/d2r;

           /* convert degrees to radians */
       x=atan(1)/45*x; 
       y=atan(1)/45*y;
           output;
    end;
  run;

DATA usotag;
  SET maps.states;
  /* PERSONALIZE: Add your states in here: If you want other countries, Add them too */ 
  IF state in (1,5,9,10,11,12,13,17,18,19,20,21,22,23,24,25,26,27,28,29,
               31,33,34,36,37,38,39,40,42,44,45,46,47,48,50,51,54,55,75);
RUN;

   
  /* combine the annotate data set with the map data set */
  data combo;
    set anno annoleg usotag;
  run; 

  /* project the data */
  proc gproject data=combo out=proj dupok;
    id state;
  run;

  /* create a new map data set and annotate data set */
  /* from the projected data */
  data tx label;
    set proj;
   country = "USA" ; 
    /* annotate data set observation */
    if flag=1 then output label;
    /* map data set observation */
    else output tx;
  run;


  DATA _NULL_;
    CALL SYMPUT ('dayname', SYSGET("EMSDAY"));
    CALL SYMPUT ('casename', SYSGET("CASE"));
    STOP;
  RUN;

pattern1 value=msolid color=white;
pattern2 value=msolid color=ulig;
pattern3 value=msolid color=paoy;
pattern4 value=msolid color=ltblue;
pattern5 value=msolid color=ltgreen;

/* PERSONALIZE: Of course, You can change the Title. Do that here */ 

title1 h=2 "Circle Plot of &polid Sources";
title2 h=1 "CASE: &casename";
FOOTNOTE j=l h=0.5 "&foot"; 

  proc gmap data=tx map=tx anno=label;
    id state; 
    choro country / coutline=black nolegend;
  run;
  quit;