*
; 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 */ *;