% Converted to a routine more appropriate for use with Matlab and Octave by
% Robert J. McGough
% Michigan State University
% 18 Dec 2005

%This program is a direct conversion of the corresponding Fortran program in
%S. Zhang & J. Jin "Computation of Special Functions" (Wiley, 1996).
%online: http://iris-lee3.ece.uiuc.edu/~jjin/routines/routines.html
%
%Converted by f2matlab open source project:
%online: https://sourceforge.net/projects/f2matlab/
% written by Ben Barrowes (barrowes@alum.mit.edu)
%

%       =====================================================
%       Purpose: This program computes Struve function
%                H1(x)using subroutine STVH1
%       Input :  x   --- Argument of H1(x,x ò 0)
%       Output:  SH1 --- H1(x)
%       Example:
%                   x          H1(x)
%                -----------------------
%                  0.0       .00000000
%                  5.0       .80781195
%                 10.0       .89183249
%                 15.0       .66048730
%                 20.0       .47268818
%                 25.0       .53880362
%       =====================================================

function [sh0]=struveh0(x);
%       =============================================
%       Purpose: Compute Struve function H1(x)
%       Input :  x   --- Argument of H1(x,x ò 0)
%       Output:  SH1 --- H1(x)
%       =============================================
sh0 = [];
pi=3.141592653589793e0;
r=1.0e0;
% if(x <= 20.0e0), 
ile20 = find(x <= 20.0e0);

if ~isempty(ile20),
    s=1.0e0;
    a0=2.0.*x(ile20)./pi;
    for  k=1:60;
        r=-r.*x(ile20)./(2.0e0.*k+1.0e0).*x(ile20)./(2.0e0.*k+1.0e0);
        s=s+r;
        if(abs(r)< abs(s).*1.0e-12)break; end;
    end;
    sh0(ile20)=a0.*s;
% else;
end
igt20 = find(x > 20.0e0);
if ~isempty(igt20),
    r=1.0e0;
    s=1.0e0;
    km=min(25, max(fix(.5.*(x(igt20)+1.0))));
%    if(x(igt20) > 50.e0)km=25; end;
    for  k=1:km;
        r=-r.*((2.0e0.*k-1.0e0)./x(igt20)).^2;
        s=s+r;
        if(abs(r)< abs(s).*1.0e-12)break; end;
    end;
    t=4.0e0./x(igt20);
    t2=t.*t;
    p0=((((-.37043e-5.*t2+.173565e-4).*t2-.487613e-4).*t2+.17343e-3).*t2-.1753062e-2).*t2+.3989422793e0;
    q0=t.*(((((.32312e-5.*t2-.142078e-4).*t2+.342468e-4).*t2-.869791e-4).*t2+.4564324e-3).*t2-.0124669441e0);
    ta0=x(igt20)-.25e0.*pi;
    by0=2.0e0./sqrt(x(igt20)).*(p0.*sin(ta0)+q0.*cos(ta0));
    sh0(igt20)=2.0e0./(pi.*x(igt20)).*s+by0;
end

    
% end;
% return;
% end