*
;
              options nosource;
/* INTEGRAL - CALCULATES INTEGRATED FUNCTION FROM A FUNCTION GIVEN
              IN PARAMETERFORM

    Written:  June 1, 1993, revised: July 7, 1993, October 12, 1993
                                     October 20, 25, 26 and 27 1994
                                     July 14, 1995, July 25, 1998

    Developed using SAS 6.07 on VM/CMS
    Runable:  SAS 6.08 - 6.12 for Windows, SAS 6.09 under AIX Unix
    Author:   Arnold Schick
    Procs:    PROC MEANS, PROC DATASETS, PROC SORT and PROC IML (SPLINE)
    Other:    SAS MACRO language
    Macros:   only INTEGRAL
    Note:     Do not use _INITIA_, _TEMP_ or _MINMAX_ or F_L as a data
              set name. Observations from input data set with missing
              values are excluded.

    Macro Call: %INTEGRAL(DATA,RESULT,X,Y,X_1,X_2,C,M,B,SCALE,NUMB,SORT);


MACRO VARIABLE      DESCRIPTION
------------------+-------------------------------------------------
 In Request:

 DATA               Name of SAS data set with input data. If this
                    parameter is missing or blank _LAST_ is used.

 RESULT             Name of SAS data set created by this macro.
                    If missing or blank stops the %macro. 5 additional
                    variables are created:
                    Y_MODIFY : modified input function by M and B
                    Y_SOURCE : result of integration, function of area
                    FACE     : function of absolute area
                    MX       : function of main point coordinate X of area
                    MY       : function of main point coordinate Y of area

 X                  Variablename of independent variable in input SAS
                    data set DATA. If missing or blank, X is used.

 Y                  Variablename of dependent variable in input SAS
                    data set DATA. If missing or blank, Y is used.

 X_1                Value of the interval start-value. If it is
                    blank or missing, the minimum X-value is used.

 X_2                Value of the interval end-value. If it is
                    blank or missing, the maximum X-value is used.

 C                  Value of integration-constant. If missing or blank,
                    zero is used.

 M                  Value of slope-factor for correction of a
                    slope in input-funtion y'. If M is zero, blank
                    or missing, a slope correction is not done.
                    If M=9999, a slope is calculated with condition:
                    y'(x_min) = y'(x_max) with center in X=0.

 B                  Value for shifting of the input-function y'
                    in Y-direction. If missing, blank or zero, no
                    shifting is done.

 SCALE              Value of scaling-factor for the output-function y
                    If missing or blank, 1.0 is used. SCALE <> 0 !

 NUMB               Value of number of interpolation-points, which
                    are based on one unit length of 1 at function-
                    polyline-length f'(x). These are calculated:
                    NUMB divided by function-polyline-length. Over
                    the interval are used: NUMB multiplicated with
                    function-polyline-length. This method considers
                    volatile functions y', their scalings and the
                    geometrical meaning of integration.
                    If NUMB is zero, blank or missing, NUMB is set
                    to 1000 divided by the maximum axis range from
                    abscissa or ordinate. On X-axis parallel y',
                    2 interpolation-points are used, if possible.
                    If NUMB is less than or equal -2, that positive
                    number are based on the full X-axis range.
                    NUMB <= -2 or NUMB >= 0 !  NUMB < OBS => OBS := NUMB

 SORT               Name of option for sorting of input data, with
                    SORT    sorts input data by independent variable X,
                            this is the default by missing or blank
                    NOSORT  input data will not be sorted.

Example:

data one;
  do x=0 to 8*atan(1) by atan(1)/4;
     y=sin(x);
     output;
   end;
run;
*%macro integral (data,result, x, y, x_1, x_2, c, m, b, scale, numb,sort);

       %integral (one ,two);
       %integral (one ,two   , x, y,   .,   .,-1,  , .,  1   ,10   ,sort);


Reference:  Schick, A. and H. Rager: 'Integration of EPR Spectra'
            in Appl. Magn. Reson. 4, 367-375 (1993)

For more information:

 Arnold Schick

 University of Marburg
 Academic Computing Center
 Hans-Meerwein-Str.
 35032 Marburg/Lahn   Germany

 Internet: schick@hrz.uni-marburg.de

 If you  find an error-condition (it is provided 'as it is')
 please let us know about this error-condition. And when you
 have good tips for better formulation with SAS, let us also know.

 Please acknowledge the authors as the provider of this integration
 method in publications that it used. Thank you very much.
  */


%macro integral (data,result,x,y,x_1,x_2,c,m,b,scale,numb,sort);
options nosource nostimer nosymbolgen nonotes nomprint;

  %if &data =  and &result=  and &x=  and &y=  %then %goto quit_4;
  %if &data    =  or &data = . %then %let data = _LAST_ ;
  %if &result  =  or &result  = . %then %goto quit_1;
  %if &x       =  or &x    = . %then %let x = x;
  %if &y       =  or &y    = . %then %let y = y;
  %if &sort    =  or &sort = . %then %let sort = sort;
%local empty;
data _NULL_ ;
    if 0 then set &data nobs=last;
    call symput('empty',last);
    stop;
run;

%if &empty < 2 %then %goto quit_3;

data _TEMP_;
  set &data;
  where ((&x is not missing) and (&y is not missing));
run;

%if &sort = sort %then %do;
     proc sort data=_TEMP_ out=_TEMP_;
       by &x;
     run;
%end;

data f_L;
    set _temp_;
    if _N_ = 1 then do; first_x=&x; first_y=&y; output; end;
    set _temp_ nobs=last point=last;
    last_x=&x; last_y=&y; output;
    stop;
run;

data _NULL_;
  set f_L;
  retain first_x first_y;
  if _N_ = 2 then do;
               delta=last_x-first_x;
               if delta=0 then d=0; else d=1;
               call symput('first_x',first_x);  call symput('first_y',first_y);
               call symput('last_x',last_x);    call symput('last_y',last_y);
               call symput('d',d);
             end;
run;

%if (&d = 0) and (&sort = sort) %then %goto quit_3;

data _NULL_;
  k=1;
  xx_1 = symget('x_1')/1.0 ;
  if xx_1 = . then do; xx_1 = &first_x; k=0; end; call symput('x_1',xx_1);
  xx_2 = symget('x_2')/1.0 ;
  if xx_2 = . then do; xx_2 = &last_x; k=0; end;  call symput('x_2',xx_2);
  call symput('k',k);
  mm = symget('m')/1.0 ;
  if mm = 9999 then mm = -(&last_y-&first_y)/(&last_x-&first_x);
  if mm = . then mm = 0 ;           call symput('m',mm);
  cc = symget('c')/1.0 ;
  if cc = . then cc = 0 ;           call symput('c',cc);
  bb = symget('b')/1.0 ;
  if bb = . then bb = 0 ;           call symput('b',bb);
  sscale = symget('scale')/1.0 ;
  if sscale = . then sscale = 1.0 ; call symput('scale',sscale);
  if sscale = 0 then nil =  1;      call symput('nil',nil);
  nnumb  = symget('numb')/1.0;
  if nnumb = . then nnumb = 0;      call symput('numb',nnumb);
run;

%if &nil = %str(1) %then %goto quit_2;
%if &k=1 %then %do;
  data _temp_;
    set _temp_ ;
    where ((&x between &x_1 and &x_2) or (&x between &x_2 and &x_1));
  run;
%end;

proc means data=_temp_ noprint min max n;
   var &x &y;
   output out=_minmax_  min=x_min y_min  max=x_max y_max n=nnn;
run;
data _NULL_;
   set _minmax_;
   retain delta_x 0 delta_y 0 ;
   delta_x = abs(x_max - x_min);
   delta_y = abs(y_max - y_min);
   call symput('nnn',nnn);
   if &numb = 0 then do;
     if delta_x > 0 then
      if delta_x >= delta_y then call symput('numb',1000/delta_x);
       else if delta_y ^= 0 then call symput('numb',1000/delta_y);
   end;
 run;

%if &numb <= -2 %then %do;
   %let n = &numb*-1;
    %goto axisbase;
%end;

%if &numb <= 0 %then %goto quit_3;

data _initia_;
  set _temp_ nobs=last end=fino;
  keep sum;
  x_p = lag(&x);
  y_p = lag(&y);
  if _N_ = 1 then sum=0;
             else sum + sqrt((&x-x_p)**2 + (&y-y_p)**2);
  if fino then do;
    if fuzz(sum) = fuzz(&last_x - &first_x) or
       fuzz(sum) = fuzz(&last_y - &first_y) or
       last = 2 then sum=1/&numb;
    if fuzz(sum) = fuzz(sqrt((&last_x-&first_x)**2 + (&last_y-&first_y)**2))
      then sum=(&nnn-1)/&numb;
    output;
  end;
run;

data _null_;
  set _initia_;
  numb = &numb;
  points = round(sum*&numb+0.5,1);
  call symput ('n',points);
  put '        note: Function-polyline-length 1 obtained ' numb 6. ;
  put '              interpolation-points, if necessary.' ;
  put '              Used points over the interval: ' points '.';
run;

%if &n < 2 %then %goto quit_3;

%axisbase : ;

proc iml;
  use _temp_;
  read all;
  call spline(yy,&x||&y,,,&n);
  xx=yy(|,1|); yy=yy(|,2|);
  create &result var{xx yy};
  append;
quit;

data &result;
* This is the core  ;
  set &result nobs=n end=last;
  length face y_source y_sourc2 y_modify Ax Ay Mx My 8 ;
  x_previo = lag(xx);
  xm = (xx-x_previo)/2 + x_previo;
  y_previo = lag(yy) + x_previo*&m + &b ;
  y_modify = yy + xx*&m + &b ;
  retain x_begin ;
  if _N_ = 1 then do;
                y_sourc2 = &c;
                y_source = 0;
                face     = 0;
                x_begin  = xx;
                Ax = 0;  Ay = 0;
                Mx = xx; My = y_modify/2;
             end;
             else do;
                y_sourc2 = (xm-x_previo)*y_previo + (xx-xm)*y_modify;
                if n >= 500 then y_sourc2 + y_sourc2 * 4*atan(1)/((n-1)**2);
                face + abs(y_sourc2);
                if (y_modify + y_previo)=0
                   then x0=xm;
                   else x0=(xx-(xx-x_previo)/3*(y_modify+2*y_previo)/(y_modify+y_previo));
                Ax + abs(y_sourc2)*x0;
                mm=(y_modify-y_previo)/2/(xx-x_previo);
                bb= y_modify/2-mm*xx;
                Ay + abs(y_sourc2)*(mm*x0+bb);
                Mx = Ax/face;
                My = Ay/face;
             end;
  y_source + y_sourc2*&scale ;
  if last then do;
    file print;
    x_end = xx ;
    m=&m; c=&c; b=&b; x_1=&x_1; x_2=&x_2; scale=&scale;
    put ' Results of Integration';
    put ' ';
    put '    area=' face ' with main point at: ' "&x" '=' Mx ' and ' "&y" '=' My;
    put '    integrated between ' x_begin ' and ' x_end ;
    put '    defined interval:  ' x_1 '  to ' x_2;
    put '    with ' n ' interpolation-points' ;
    put '    input-function corrected with slope-factor=' m
        ' and shifted by ' b ;
    put '    integration constant=' c ' and scaling-factor='
          scale ' of output-function.';
  end;
  keep xx yy y_modify y_source face Mx My;
run;

proc datasets nolist;
  delete _initia_ _minmax_ _temp_ f_L;
quit;

options notes;
  data &result;
    set &result;
    rename xx=&x yy=&y;
  run;
options nonotes;

%goto final ;
%quit_1 : %put MACRO-HALT: Please define result data set. ;
  %goto final;
%quit_2 : %put MACRO-HALT: Scaling-factor equal zero. ;
  proc datasets nolist;
    delete _temp_ f_L;
  quit;
  %goto final;
%quit_3 : %put MACRO-HALT: Interpolation-points less than 2
 or values not in;
          %put defined interval or interval
 equal zero or input data set is empty. ;
%goto final;
%quit_4 : %put MACRO-HALT: Please call with possible parameters:
 integral(data,result,x,y,lower,upper,constant,slope,shift,scale,number,sort);
%final : ;
%put;
options source stimer notes;
%mend integral;   options source;


/* Example
data one;
  do x=0 to 4*atan(1) by atan(1)/8;
     y=sin(x);
     output;
  end;
run;

*%macro integral (data,result,x,y,x_1,x_2, c , m  , b  , scale,numb,sort);
  %integral      (one ,two   , , ,   ,   , 0 ,    ,    ,      ,-200,nosort);
*/


*Utiliy Macro to remark test output for example in macro INTEGRAL above;

       options nosource;
%macro remarks (data,result,x,y,text,style,size,position);
       options nostimer nonotes nomprint nosymbolgen;

*written: June 26, 1995, Arnold Schick University of Marburg;

  %if &data =  and &result =  and &x =  and &y =  or
     %upcase(&data) = HELP or &data = ? %then %do;
     %put  This macro creates annotate data set;
     %put  for annotation of text within graphs;
     %put  and prepares axes-statements: horizontal: axis1, vertical: axis2;
     %put;
     %put  call: remarks (data,result,x,y,text,font,fontsize,position);
     %put;
     %put  with meanings and defaults: ;
     %put  data     : input data set [ _last_ ];
     %put  result   : output data set [ _anno_ ];
     %put  x        : variablename for abscissa in input data set [ x ];
     %put  y        : variablename for ordinate in input data set [ y ];
     %put  text     : variablename of textvariable [ text ] in data set;
     %put  font     : variablename of font [ style ] or fontname [ simplex ];
     %put  fontsize : variablename of fontsize [ size ] or fontsize [ 1 ];
     %put  position : variablename of position [ position ] or position [ 2 ] = center;
     %put  Please note, text should be defined in format $char30;
     %put ;
     %goto fin;
  %end;
  %if &data   =. %then %let data   = _LAST_ ;
  %if &data   =  %then %goto quit_1;
  %if &result = or &result = . %then %let result = _ANNO_  ;
  %if &x =  or &x = . %then %let x = x;
  %if &y =  or &y = . %then %let y = y;
  %if &text     =  or &text=.  %then %let text =text;
  %if &style    =  or &style=. %then %let style=simplex;
  %if &size     =  %then %let size= . ;
  %if &position =  or &position =. %then %let position=2;

data &result;
    set &data;
    retain function 'LABEL' xsys '2' ysys '2' style "&style"
           size &size position "&position" text ' ';
    keep function x y xsys ysys style size position text;
    x = &x;
    y = &y;
    text = &text;
run;

proc means data=&data noprint;
  var &x &y;
  output out=_minmax_ min=min_x min_y max=max_x max_y;
run;
%let prop = %scan(&size,1,1234567890);
%if &prop ^=. %then %do;
data _NULL_;
  set _MINMAX_;
  dx = (max_x-min_x)*0.1;
  dy = (max_y-min_y)*0.1;
  call symput('min_x',min_x-dx*&size); call symput('min_y',min_y-dy*&size);
  call symput('max_x',max_x+dx*&size); call symput('max_y',max_y+dy*&size);
  call symput('dx',dx);                call symput('dy',dy);
run;
%end; %else %do;
data _NULL_;
  set _MINMAX_;
  dx = (max_x-min_x)*0.1;
  dy = (max_y-min_y)*0.1;
  call symput('min_x',min_x-dx); call symput('min_y',min_y-dy);
  call symput('max_x',max_x+dx); call symput('max_y',max_y+dy);
  call symput('dx',dx);          call symput('dy',dy);
run;
%end;

proc datasets nolist;
  delete _MINMAX_;
quit;

axis1 order=&min_x to &max_x by &dx;
axis2 order=&min_y to &max_y by &dy;

%goto fin;
%quit_1 : MACRO-HALT: %put input data set name is not present;
%fin : ;

options source stimer ;
%mend remarks;  options source;

/*
*Example for macro REMARKS;

data one;
  input x y size position $1. style $8. text $char30.;
cards;
1 1 1.1 3 zapf    Text 1
2 2 0.6 2 simplex noch ein Text
3 2 1.2 4 zapfb   ein weiterer Text
3 3   2 7 brush   letzter Text
;
run;

%remarks(one,remarks);

proc gplot data=one;
  symbol i=none v=star r=123;
  plot y*x / annotate=remarks
             haxis=axis1 vaxis=axis2
             nolegend noaxis;
run; quit;
*END of Example for macro REMARKS;
*/

/* Test Example for macro INTEGRAL;  */
data one;
  do x=0 to 4*atan(1) by atan(1)/8;
     y=sin(x);
     output;
  end;
run;

*%macro integral (data,result,x,y,x_1,x_2, c , m  , b  , scale,numb,sort);
  %integral      (one ,two   , , ,   ,   , 0 ,    ,    ,      ,-200,nosort);

data comments;
  input x y size position $1. style $8. text $char30.;
cards;
3   0.1 1.0 3 zapf    Original Function f(x)
3.2 1.8 1.0 2 zapf    Integral of f(x)
;
run;

 %remarks(comments,remarks);

proc gplot data=two;
  symbol i=join r=2;
  plot (y y_source) * x / overlay annotate=remarks;
run; quit;
 /* End of example */
 *
;