function y=powerdist(x,alpha,ymin,ymax)
% convert array of uniform deviate random numbers to power law distribution 
% we aim to use fomulae:
%    F(x) = ( (ymax^(1-alpha) - ymin^(1-alpha) ) * x + ymin^(1-alpha) )^(1./(1-alpha))
%           
%
% x     - array of unirom random numbes in the range of [0,1]
% ymin    - cutoff value for power law distribution
% ymax  - maximum value for power law distribution 
% alpha - power law exponent
% F(x)  - array of powerlaw distribution with exponent -alpha 
% safari@znu.ac.ir
%% when you publish results based on this code, please refer to the 
%% following paper: 
% %% 
% @ARTICLE{2012ApJ...744..113T,
%    author = {{Tajfirouze}, E. and {Safari}, H.},
%     title = "{Can a Nanoflare Model of Extreme-ultraviolet Irradiances Describe the Heating of the Solar Corona?}",
%   journal = {\apj},
% archivePrefix = "arXiv",
%    eprint = {1112.4926},
%  primaryClass = "astro-ph.SR",
%  keywords = {Sun: corona, Sun: flares, Sun: UV radiation},
%      year = 2012,
%     month = jan,
%    volume = 744,
%       eid = {113},
%     pages = {113},
%       doi = {10.1088/0004-637X/744/2/113},
%    adsurl = {http://adsabs.harvard.edu/abs/2012ApJ...744..113T},
%   adsnote = {Provided by the SAO/NASA Astrophysics Data System}
% }
 
% Example:
% x=random('unif',0,1,10000,1);
% figure(1),subplot(1,2,1),plot(x(1:100))
% figure(1),subplot(1,2,2),plot(hist(x,100))
%
% alpha=2.0; ymin=0.003; ymax=0.01; 
%
% y=powerdist(x,alpha,ymin,ymax)
% 
% figure(2),subplot(1,2,1),plot(y(1:100))
% figure(1),subplot(1,2,2),plot(hist(y,100))
%  
% test for power law fit 
% [hy,yy]=hist(y,100);
%  cf=powerlaw_fit(yy,hy)
%
par1=1-alpha; 
par2=1./par1;
par3= ymax^par1-ymin^par1;
y=( par3*x+ymin^par1).^par2; 
end
