Re: Limb Darkening

Hi all,
 looks like we are heading to a nice 2014 with GMT 5 out and good discussions. A few months ago I was trying to do "scientific illustration" with an Earth with night and daysides, so sort of what is asked for here. See my script fragment below and note some symbols are undefined in this fragment, but I think you can figure them out.

 In attachment is the code for illumination. As noted in Eduardo Suarez' posting, it is hard to find a "nice decay rate". For testing, my code can be uncommented to have a rather sharp day/night boundary. The illumination is then overlaid with transparency, and to see it you have to make a PDF file.

Enjoy...Martin

paltb=gray
makecpt -C$paltb -T0/1/0.01 > illum.cpt

# center of projection
lonp=280.34
lonp=320
latp=75

projection="JG$lonp/$latp/$D1"
region="Rg"

./illumination $geonoon | surface -Rg -I1/0.1 -T0 -Gillum.grd

pscoast -$region -Dc -A10000 -$projection -B30g30/15g15 -W1 -X$xoff -Y$yoff -Sblue -Ggreen $orient -K > illum.ps

# grdimage with transparency -t% from 0 (opaque) to 100 (transparent)
trans=10
grdimage illum.grd -$region -Cillum.cpt -$projection -K -O -t$trans >> illum.ps

ps2pdf illum.ps

-- 
    This communication is intended for the use of the recipient to whom it
    is addressed, and may contain confidential, personal, and or privileged
    information. Please contact us immediately if you are not the intended
    recipient of this communication, and do not copy, distribute, or take
    action relying on it. Any communications received in error, or
    subsequent reply, should be deleted or destroyed.
---

To unsubscribe, send the message "signoff gmt-help" to [email protected]

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

#define DEGREE M_PI/180.

main(int argc, char **argv)
{
/* calculate illumination as a function of great circle distance to subsolar point entered as argmument lonss, latss */
/* this is mainly oriented to global gridding for making an overlay of illumination for scientific illustration      */
double lonss,latss,lonssr,latssr,angdistssrn,lon,lat,lonr,latr,illum;
double angdistssr(double lonssr,double latssr,double lonr,double latr);
if(argc!=3) { printf("usage %s lon_subsolar lat_subsolar\n",argv[0]); exit(0); }
lonss=atof(argv[1]); if(lonss<0.||lonss>360.) {printf("lon_subsolar (%s) should be in range [0.,360.]\n",argv[1]); exit(0); }
latss=atof(argv[2]); if(latss<-90.||latss>90.) {printf("lat_subsolar (%s) should be in range [-90.,90.]\n",argv[2]); exit(0); }
lonssr=lonss*DEGREE; latssr=latss*DEGREE;
for(lat=-90.;lat<=90;lat+=2.)
 for(lon=0.;lon<=360;lon+=15.)
  { lonr=lon*DEGREE; latr=lat*DEGREE;
    angdistssrn=angdistssr(lonssr,latssr,lonr,latr);
//    illum=(angdistssrn>=M_PI/2)? 0.0 : 1.0;
/* "logistic curve smoothed illumination near terminator */
    illum=1.0-1.0/(1+exp(-10.0*(angdistssrn-(M_PI/2.0))));
    printf("%lf %lf %lf\n",lon,lat,illum);
  }
}

double angdistssr(double lonssr,double latssr,double lonr,double latr)
{
/* from Wikipedia article on Great Circle Distance - note there can be numerical problems with close points */
return acos(sin(latr)*sin(latssr)+cos(latr)*cos(latssr)*cos(fabs(lonr-lonssr)));
}