Initial upload of filesfrom/noc/mpoc/rpdmoc/jc145/software_jc145/mexec_v3

pstar/progs/Contents.m

+% Selection of programs to manipulate pstar binary files in matlab
+% 
+% plisth  : load a pstar file header and echo it to the screen
+% pload   : load a pstar file; fixes variable names that would be illegal in matlab. (calls plisth).
+% pst2nc  : copy contents of pstar file to netcdf file (calls pload).
+% nclisth : list the header of a pstar to netcdf file made by pst2nc
+%

pstar/progs/plisth.m

+function h = plisth(filename)
+
+% plisth: load a pstar file header and echo it to the screen
+%use:
+% h = plisth(filename)
+%where filename is the name of a pstar file.
+%
+% BAK at SOC 31 March 2005
+
+if nargin < 1
+    filename = input('Type name of pstar file to load :\n','s');
+end
+
+if exist(filename,'file') ~= 2
+    disp(['Filename ''' filename ''' not found']);
+    return
+end
+
+magvar=1344950870;
+magdat=1344950852;
+magvr8=1344951864;
+magdr8=1344947256;
+
+r5 = 1365;
+r8 = 1024;
+
+fid = fopen(filename,'r','b');
+
+h.datnam = char(fread(fid,8,'*uchar'))';
+h.vers = char(fread(fid,2,'*uchar'))';
+h.opwrit = char(fread(fid,1,'*uchar'))';
+h.rawdat = char(fread(fid,1,'*uchar'))';
+h.pipefl = char(fread(fid,1,'*uchar'))';
+h.archiv = char(fread(fid,1,'*uchar'))';
+h.magic = fread(fid,1,'uint32');
+
+lb = 0;
+if h.magic == magvar | h.magic == magdat 
+    lb = r5;
+elseif h.magic == magvr8 | h.magic == magdr8
+    lb = r8;
+else
+    disp('Problem with magic number - quitting')
+    return
+end
+
+storvar = '   ';
+if h.magic == magvar | h.magic == magvr8
+    storvar = 'var';
+else
+    storvar = 'dcs';
+end
+
+blank6 = char(fread(fid,6,'*uchar'))';
+h.prefil = char(fread(fid,8,'*uchar'))';
+h.postfl = char(fread(fid,8,'*uchar'))';
+h.noflds = fread(fid,1,'uint32');
+h.norecs = fread(fid,1,'uint32');
+h.nrows = fread(fid,1,'uint32');
+h.nplane = fread(fid,1,'uint32');
+h.icent = fread(fid,1,'uint32');
+h.iymd = fread(fid,1,'uint32');
+h.ihms = fread(fid,1,'uint32');
+h.platnam = char(fread(fid,12,'*uchar'))';
+h.platyp = char(fread(fid,8,'*uchar'))';
+h.pltnum = char(fread(fid,8,'*uchar'))';
+h.instmt = char(fread(fid,12,'*uchar'))';
+blank4 = char(fread(fid,4,'*uchar'))';
+h.recint = char(fread(fid,16,'*uchar'))';
+
+if lb == r5 %read them as real*5 and unpack; skip the status byte at read time. Then skip 2 blanks.
+    alat6 = char(fread(fid,5,'5*uchar',1))';
+    blank2 = char(fread(fid,2,'*uchar'))';
+    alat = punpack(alat6);
+    along6 = char(fread(fid,5,'5*uchar',1))';    
+    blank2 = char(fread(fid,2,'*uchar'))';
+%     blank2 = char(fread(fid,2,'5*uchar',1))';
+    % long standing bug in above line fixed by BAK on 26 Aug 2008.
+    % In the above line, there is an incorrect (not required) skip as the fourth argument of the call to fread.
+    % In older versions of matlab, the fact that 2 bytes were read to 5*uchar meant that the skip was
+    % never executed. In a newer version running on solaris workstation 'rapid', the skip was executed once
+    % even though the 5*uchar was not filled. This meant an extra byte was read out of the pstar header,
+    % and all following data were corrupt, especially the filenames. This error was obviously apparent,
+    % far too obvious to inadvertantly not notice. So there is no danger of loaded files 'aprearing to be OK but actually bad'
+    along = punpack(along6);
+    dpthi6 = char(fread(fid,5,'5*uchar',1))';
+    blank2 = char(fread(fid,2,'*uchar'))';
+    dpthi = punpack(dpthi6);
+    dpthw6 = char(fread(fid,5,'5*uchar',1))';
+    blank2 = char(fread(fid,2,'*uchar'))';
+    dpthw = punpack(dpthw6);
+else %lb == r8; read them as real*8
+    alat = fread(fid,1,'double');
+    along = fread(fid,1,'double');
+    dpthi = fread(fid,1,'double');
+    dpthw = fread(fid,1,'double');
+end
+h.alat = alat;
+h.along = along;
+h.dpthi = dpthi;
+h.dpthw = dpthw;
+
+for k = 1:12
+    h.coment(k,1:72) = char(fread(fid,72,'*uchar'))';
+end
+ifldxx = 128;
+for k = 1:ifldxx
+    h.fldnam{k} = char(fread(fid,8,'*uchar'))';
+end
+
+for k = 1:ifldxx
+    h.fldunt{k} = char(fread(fid,8,'*uchar'))';
+end
+
+if lb == r5
+    %128*5 = 640; We need 640 bytes in fives, which gives 128 5-byte elements.
+    %we skip 1 status byte with each element, so we read 128*(5+1) = 768 bytes.
+    %So skip 256 more, so we advance 1024 bytes for each array.
+    alrlim_array = char(fread(fid,ifldxx*5,'5*uchar',1))';
+    fseek(fid,256,0);
+    uprlim_array = char(fread(fid,ifldxx*5,'5*uchar',1))';
+    fseek(fid,256,0);
+    absent_array = char(fread(fid,ifldxx*5,'5*uchar',1))';
+    fseek(fid,256,0);
+    alrlim = punpack(alrlim_array(1:ifldxx*5));
+    uprlim = punpack(uprlim_array(1:ifldxx*5));
+    absent = punpack(absent_array(1:ifldxx*5));
+else %lb == r8
+    alrlim = fread(fid,128,'double');
+    uprlim = fread(fid,128,'double');
+    absent = fread(fid,128,'double');
+end
+
+h.alrlim = alrlim(:)';
+h.uprlim = uprlim(:)';
+h.absent = absent(:)';
+
+blank = char(fread(fid,2046,'*uchar'))';
+h.site = char(fread(fid,2,'*uchar'))';
+
+h
+
+latd = fix(alat);
+latm = abs(60*(alat-latd));
+lond = fix(along);
+lonm = abs(60*(along-lond));
+disp(['Data Name :  ' h.datnam ' ' h.site ' ' h.vers]);
+disp(['Platform :   ' h.platyp ' ' h.platnam ' ' h.pltnum]);
+disp(['Instrument :' h.instmt '   dpthi ' sprintf('%8.2f',h.dpthi) '   dpthw ' sprintf('%8.2f',h.dpthw)]);
+disp(['Fields :    ' sprintf('%3d',h.noflds) '     Data Cycles ' sprintf('%8d',h.norecs) '     Rows ' sprintf('%4d',h.nrows) '   Planes ' sprintf('%4d',h.nplane)]);
+disp(['Position (lat lon) : '  sprintf('%10.5f',alat) '  ' sprintf('%10.5f',along)]);
+disp(['Position (lat lon) : '  sprintf('%4d %06.3f',latd,latm) ' ' sprintf('%4d %06.3f',lond,lonm)]);
+disp(['Time origin : ' sprintf('%02d',h.icent/100) '/' sprintf('%06d',h.iymd) '/' sprintf('%06d',h.ihms)]);
+disp('************************************************************************');
+for k = 1:h.noflds
+    disp(['*' sprintf('%3d',k) '*' h.fldnam{k} '*' h.fldunt{k} '* ' sprintf('%15.3f',h.alrlim(k)) ' * ' sprintf('%15.3f',h.uprlim(k)) ' * ' sprintf('%10.3f',h.absent(k)) ' *']);
+end
+disp('************************************************************************');
+for k = 1:12
+    thiscoment = h.coment(k,:);
+    blank72 = '                                                                        ';
+    if strcmp(blank72,thiscoment) == 0
+        disp(thiscoment);
+    end
+end
+
+
+fclose(fid);
+
+return

pstar/progs/pload.m

+function [d, h] = pload(filename,instring)
+
+% pload: load a pstar file; fixes variable names that would be illegal in matlab.
+% filename is a pstar filename
+% instring is an optional character array of variable numbers/names to load; use '/' or '' or ' ' for  all;
+% [d h] = pload(filename,instring);   %full form of command
+% [d h] = pload;              %will prompt for filename and varlist
+% [d h] = pload(filename);    %will prompt for variable list
+% [d h] = pload(filename,''); %will select all variables
+% [d h] = pload(filename,'1 3~5 a b');  %will select vars 1 3,4,5 and variables named 'a' and 'b'
+%
+% The final list of variable numbers is passed through 'unique' before loading.
+%
+% the variables are loaded into struct array d
+% the header info is loaded into struct array h
+%
+% BAK at SOC 31 March 2005
+
+% As presently written, the load is a bit slow for pstar data packed as
+% real*5; The time is mostly spent in punpack.m
+
+
+clear d h
+
+if nargin < 1
+    filename = input('Type name of pstar file to load :\n','s');
+elseif nargin > 2
+    disp(['You must give a maximum of two arguments']);
+    return
+end
+
+if exist(filename,'file') ~= 2
+    disp(['Filename ''' filename ''' not found']);
+    return
+end
+
+h = plisth(filename);
+norecs = h.norecs;
+noflds = h.noflds;
+nrows = max(h.nrows,1);
+%bak at NOCS, 14 Jul 2005: include possibility of NPLANE; I didn't know anyone actually used this,
+% but RTP says he sometimes does !
+nplane = max(h.nplane,1);
+
+if mod(norecs,nrows*nplane) ~= 0
+    error(['nrows * nplane ' num2str(nrows) ' * ' num2str(nplane) ' does not divide norecs ' num2str(norecs)])
+    return
+else
+    ncols = round(norecs/(nrows*nplane));
+end 
+
+if nargin < 2
+    disp([' ']);
+    disp(['Type list of pstar variables to load']);
+    disp(['List can be space or comma separated']);
+    instring = input('Use variable names or numbers or n1~n2. Type c/r or / for all); \n','s');
+end
+
+vlist = getvlist(instring,h);
+
+vlist = unique(vlist); %sort and remove duplicates for load
+
+skip = 1+0*[1:noflds];
+skip(vlist) = 0;
+loadlist = vlist;
+
+magvar=1344950870;
+magdat=1344950852;
+magvr8=1344951864;
+magdr8=1344947256;
+
+r5 = 1365;
+r8 = 1024;
+
+lb = 0;
+if h.magic == magvar | h.magic == magdat 
+    lb = r5;
+elseif h.magic == magvr8 | h.magic == magdr8
+    lb = r8;
+else
+    disp('Problem with magic number - quitting')
+    return
+end
+
+storvar = '   ';
+if h.magic == magvar | h.magic == magvr8
+    storvar = 'var';
+else
+    storvar = 'dcs';
+end
+
+varnames = '';
+disp(' ');
+printsummary = 0;  %Will set this to 1 if any variables are renamed
+for k = 1:noflds %sort out difficult variable names
+    if skip(k) == 1
+        varnames{k} = '';
+        continue
+    end
+    vn = h.fldnam{k};
+    vnold = vn;
+    renamed = 0;
+    %list of problem characters that may occur in pstar file names
+    %you can add to the list as new ones are discovered.
+    %Blanks are simply discarded.
+    if strcmp(vn,'        ') % if the pstar fldnam was all blanks, rename the variable.
+        vn = ['blank_' num2str(k) '_'];
+        renamed = 1;
+    end
+    
+    %remove trailing blanks;
+    while strcmp(vn(end),' ')
+        vn(end) = [];
+    end
+    
+    swap = {
+        ' ' '_space_'
+        '-' '_minus_'
+        '+' '_plus_'
+        '/' '_slash_'
+        '*' '_star_'
+        '.' '_dot_'
+        '#' '_hash_'
+        '$' '_dollar_'
+        '^' '_hat_'
+        '&' '_amp_'
+        '(' '_lparen_'
+        ')' '_rparen_'
+        '[' '_lbrac_'
+        ']' '_rbrac_'
+    } ;
+    nswaps = size(swap,1);
+    for k2 = 1:nswaps
+        s1 = char(swap(k2,1));
+        s2 = char(swap(k2,2));
+        sindex = findstr(vn,s1);
+        if length(sindex) > 0
+            vnew = '';
+            count = 1;
+            for kswap = 1:length(sindex)
+                vnew = [vnew vn(count:sindex(kswap)-1) s2];
+                count = sindex(kswap)+1;
+            end
+            vnew = [vnew vn(count:end)];
+            vn = vnew;
+            renamed = 1;
+        end
+    end
+    
+    %replace any other odd characters with underscore. We'll permit alphanumeric and underscore in matlab var names.
+    okchars = '_0123456789abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ';
+    for k2 = 1:length(vn)
+        if length(findstr(okchars,vn(k2))) == 0
+            vn(k2) = '_';
+            renamed = 1;
+        end
+    end
+
+    okfirstchar = 'abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ'; %must start with an alpha; if not, insert 'v_'
+    if length(findstr(okfirstchar,vn(1))) == 0
+        vn = ['v_' vn];
+        renamed = 1;
+    end
+    
+    if length(vn) > 26 %truncate. Allow for later addition of _nnn_ and still keep
+                       %total var name to shorter than 31, which seems to be needed in 
+                       % Matlab fieldnames in a structure.
+        vn = vn(1:26);
+        renamed = 1;
+    end
+    
+    
+    smatch = strmatch(vn,varnames,'exact');
+    if length(smatch) > 0
+        vn = [vn '_' num2str(k) '_'];
+        renamed = 1;
+    end
+    
+    if renamed == 1
+        printsummary = 1;
+        disp(['Variable number ' sprintf('%3d',k) ':''' vnold ''' renamed to ''' vn '''']);     
+    end
+    
+    vnprint = vn;
+    while length(vnprint) < 32
+        vnprint = [vnprint ' '];
+    end
+    varnames{k} = vnprint;
+end % of for k = 1:noflds %sort out problem variable names
+disp(' ');
+if printsummary == 1
+    disp('List of loaded variables will be');
+    disp('************************************************************************************************');
+    for k = loadlist
+        disp(['*' sprintf('%3d',k) '*' sprintf('%32s',varnames{k}) '*' h.fldunt{k} '* ' sprintf('%15.3f',h.alrlim(k)) ' * ' sprintf('%15.3f',h.uprlim(k)) ' * ' sprintf('%10.3f',h.absent(k)) ' *']);
+    end
+    disp('************************************************************************************************');
+    disp(' ')
+end
+
+if strcmp(storvar,'var') %store by variable
+    numblocks_per_var = ceil(norecs/lb);
+    num_data_blocks = noflds*numblocks_per_var;
+else %store by data cycle
+    numcycles_per_block = floor(lb/noflds);
+    num_data_blocks = ceil(norecs/numcycles_per_block);
+end
+
+disp(['There are ' num2str(num_data_blocks) ' data blocks to read'])
+
+fid = fopen(filename,'r','b');
+fseek(fid,8192,0); %skip header block
+
+% reserve space for loaded variables to speed execution in loop
+for kv = loadlist
+    vn = varnames{kv};
+    vn = m_remove_outside_spaces(vn);
+    cmd = ['d.' vn ' = nan+ones(1,h.norecs);']; eval(cmd)
+end
+
+for k = 1:num_data_blocks
+    if mod(k,100) == 0;
+        disp(['Reading data block ' num2str(k) ' out of ' num2str(num_data_blocks)]);
+    end;
+
+    clear buf raw
+    numvar = 1 + mod(k-1,noflds); %this is the numvar if storage is by variable
+    if strcmp(storvar,'var') & skip(numvar) == 1
+        fseek(fid,8192,0); %skip a block
+        continue
+    elseif lb == r8;
+        buf = fread(fid,lb,'double'); %read from real*8
+    else
+        %1365*5 = 6825
+        raw = char(fread(fid,6825,'5*uchar',1)); %read some real*5; skip the status byte
+        % The syntax of fread means you collect 6825 bytes of data in 5s, (ie 1365 elements) and 1365 skips of 1 byte.
+        % So there are 2 bytes left in the block of 8192.
+        fseek(fid,2,0); %skip the last 2 bytes
+        buf = punpack(raw(1:6825));
+    end
+
+    %now store the data
+
+    if strcmp(storvar,'var') %store this block in the correct variable
+        blocknum_thisvar = 1+floor((k-1)/noflds);
+        completed = lb*(blocknum_thisvar-1);
+        remaining = norecs-completed;
+        thistime = min(lb,remaining);
+        dcnums = completed+1:completed+thistime;
+        vn = varnames{numvar};
+        eval(['d.' vn '(dcnums) = buf(1:thistime);']);
+    else %store by data cycle; distribute data across variables
+        completed = numcycles_per_block*(k-1);
+        remaining = norecs-completed;
+        for k2 = 1:min(numcycles_per_block,remaining)
+            for k3 = loadlist
+                bufindex = noflds*(k2-1)+k3;
+                vn = varnames{k3};
+                eval(['d.' vn '(completed+k2) = buf(bufindex);']);
+            end %of k3 = 1:noflds
+        end %of k2 = 1:numcycles_per_block
+    end %of assigning data to variables
+
+end %of k = 1:num_data_blocks
+disp(['Finished reading data blocks']);
+
+fclose(fid);
+
+for k3 = loadlist  %sort out gridded data and absent values
+    clear grid vec absent bad
+    vn = varnames{k3};
+    eval(['vec = d.' vn ';']);
+    absent = h.absent(k3);
+    bad = find(vec == absent);
+    vec(bad) = nan;  %set absent data to nan;
+    if nrows > 1 & nplane < 2
+        vec = reshape(vec,nrows,ncols);
+    end
+    if nplane > 1
+        vec = reshape(vec,nrows,ncols,nplane);
+    end
+    eval(['d.' vn ' = vec;']);
+    %if nrows > 1 %assign data to 2-D matrix %better done with reshape command as above
+    %    for k1 = 1:nrows
+    %        for k2 = 1:ncols
+    %            grid(k1,k2) = vec(nrows*(k2-1)+k1);
+    %        end
+    %    end
+    %    eval(['d.' vn ' = grid;']); 
+    %else %copy one-D matrix with NaNs back to d.vn
+    %    eval(['d.' vn ' = vec;']);
+    %end
+end
+
+function sot = m_remove_outside_spaces(sin)
+
+% remove any outside spaces from a string
+
+s = [ ' ' sin ' '];
+
+while strcmp(s(1),' ') == 1
+    s(1) = [];
+    if isempty(s); break; end
+end
+
+s = [s ' '];
+while strcmp(s(end),' ') == 1
+    s(end) = [];
+    if isempty(s); break; end
+end
+
+sot = s;
+return
+

pstar/progs/pload_old.m

+function [d, h] = pload(filename,instring)
+
+% pload: load a pstar file; fixes variable names that would be illegal in matlab.
+% filename is a pstar filename
+% instring is an optional character array of variable numbers/names to load; use '/' or '' or ' ' for  all;
+% [d h] = pload(filename,instring);   %full form of command
+% [d h] = pload;              %will prompt for filename and varlist
+% [d h] = pload(filename);    %will prompt for variable list
+% [d h] = pload(filename,''); %will select all variables
+% [d h] = pload(filename,'1 3~5 a b');  %will select vars 1 3,4,5 and variables named 'a' and 'b'
+%
+% The final list of variable numbers is passed through 'unique' before loading.
+%
+% the variables are loaded into struct array d
+% the header info is loaded into struct array h
+%
+% BAK at SOC 31 March 2005
+
+% As presently written, the load is a bit slow for pstar data packed as
+% real*5; The time is mostly spent in punpack.m
+
+
+clear d h
+
+if nargin < 1
+    filename = input('Type name of pstar file to load :\n','s');
+elseif nargin > 2
+    disp(['You must give a maximum of two arguments']);
+    return
+end
+
+if exist(filename,'file') ~= 2
+    disp(['Filename ''' filename ''' not found']);
+    return
+end
+
+h = plisth(filename);
+norecs = h.norecs;
+noflds = h.noflds;
+nrows = max(h.nrows,1);
+%bak at NOCS, 14 Jul 2005: include possibility of NPLANE; I didn't know anyone actually used this,
+% but RTP says he sometimes does !
+nplane = max(h.nplane,1);
+
+if mod(norecs,nrows*nplane) ~= 0
+    error(['nrows * nplane ' num2str(nrows) ' * ' num2str(nplane) ' does not divide norecs ' num2str(norecs)])
+    return
+else
+    ncols = round(norecs/(nrows*nplane));
+end 
+
+if nargin < 2
+    disp([' ']);
+    disp(['Type list of pstar variables to load']);
+    disp(['List can be space or comma separated']);
+    instring = input('Use variable names or numbers or n1~n2. Type c/r or / for all); \n','s');
+end
+
+vlist = getvlist(instring,h);
+
+vlist = unique(vlist); %sort and remove duplicates for load
+
+skip = 1+0*[1:noflds];
+skip(vlist) = 0;
+loadlist = vlist;
+
+magvar=1344950870;
+magdat=1344950852;
+magvr8=1344951864;
+magdr8=1344947256;
+
+r5 = 1365;
+r8 = 1024;
+
+lb = 0;
+if h.magic == magvar | h.magic == magdat 
+    lb = r5;
+elseif h.magic == magvr8 | h.magic == magdr8
+    lb = r8;
+else
+    disp('Problem with magic number - quitting')
+    return
+end
+
+storvar = '   ';
+if h.magic == magvar | h.magic == magvr8
+    storvar = 'var';
+else
+    storvar = 'dcs';
+end
+
+varnames = '';
+disp(' ');
+printsummary = 0;  %Will set this to 1 if any variables are renamed
+for k = 1:noflds %sort out difficult variable names
+    if skip(k) == 1
+        varnames{k} = '';
+        continue
+    end
+    vn = h.fldnam{k};
+    vnold = vn;
+    renamed = 0;
+    %list of problem characters that may occur in pstar file names
+    %you can add to the list as new ones are discovered.
+    %Blanks are simply discarded.
+    if strcmp(vn,'        ') % if the pstar fldnam was all blanks, rename the variable.
+        vn = ['blank_' num2str(k) '_'];
+        renamed = 1;
+    end
+    
+    %remove trailing blanks;
+    while strcmp(vn(end),' ')
+        vn(end) = [];
+    end
+    
+    swap = {
+        ' ' '_space_'
+        '-' '_minus_'
+        '+' '_plus_'
+        '/' '_slash_'
+        '*' '_star_'
+        '.' '_dot_'
+        '#' '_hash_'
+        '$' '_dollar_'
+        '^' '_hat_'
+        '&' '_amp_'
+        '(' '_lparen_'
+        ')' '_rparen_'
+        '[' '_lbrac_'
+        ']' '_rbrac_'
+    } ;
+    nswaps = size(swap,1);
+    for k2 = 1:nswaps
+        s1 = char(swap(k2,1));
+        s2 = char(swap(k2,2));
+        sindex = findstr(vn,s1);
+        if length(sindex) > 0
+            vnew = '';
+            count = 1;
+            for kswap = 1:length(sindex)
+                vnew = [vnew vn(count:sindex(kswap)-1) s2];
+                count = sindex(kswap)+1;
+            end
+            vnew = [vnew vn(count:end)];
+            vn = vnew;
+            renamed = 1;
+        end
+    end
+    
+    %replace any other odd characters with underscore. We'll permit alphanumeric and underscore in matlab var names.
+    okchars = '_0123456789abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ';
+    for k2 = 1:length(vn)
+        if length(findstr(okchars,vn(k2))) == 0
+            vn(k2) = '_';
+            renamed = 1;
+        end
+    end
+
+    okfirstchar = 'abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ'; %must start with an alpha; if not, insert 'v_'
+    if length(findstr(okfirstchar,vn(1))) == 0
+        vn = ['v_' vn];
+        renamed = 1;
+    end
+    
+    if length(vn) > 26 %truncate. Allow for later addition of _nnn_ and still keep
+                       %total var name to shorter than 31, which seems to be needed in 
+                       % Matlab fieldnames in a structure.
+        vn = vn(1:26);
+        renamed = 1;
+    end
+    
+    
+    smatch = strmatch(vn,varnames,'exact');
+    if length(smatch) > 0
+        vn = [vn '_' num2str(k) '_'];
+        renamed = 1;
+    end
+    
+    if renamed == 1
+        printsummary = 1;
+        disp(['Variable number ' sprintf('%3d',k) ':''' vnold ''' renamed to ''' vn '''']);     
+    end
+    
+    vnprint = vn;
+    while length(vnprint) < 32
+        vnprint = [vnprint ' '];
+    end
+    varnames{k} = vnprint;
+end % of for k = 1:noflds %sort out problem variable names
+disp(' ');
+if printsummary == 1
+    disp('List of loaded variables will be');
+    disp('************************************************************************************************');
+    for k = loadlist
+        disp(['*' sprintf('%3d',k) '*' sprintf('%32s',varnames{k}) '*' h.fldunt{k} '* ' sprintf('%15.3f',h.alrlim(k)) ' * ' sprintf('%15.3f',h.uprlim(k)) ' * ' sprintf('%10.3f',h.absent(k)) ' *']);
+    end
+    disp('************************************************************************************************');
+    disp(' ')
+end
+
+if strcmp(storvar,'var') %store by variable
+    numblocks_per_var = ceil(norecs/lb);
+    num_data_blocks = noflds*numblocks_per_var;
+else %store by data cycle
+    numcycles_per_block = floor(lb/noflds);
+    num_data_blocks = ceil(norecs/numcycles_per_block);
+end
+
+disp(['There are ' num2str(num_data_blocks) ' data blocks to read'])
+
+fid = fopen(filename,'r','b');
+fseek(fid,8192,0); %skip header block
+
+for k = 1:num_data_blocks
+    if mod(k,10) == 0; disp(['Reading data block ' num2str(k) ' out of ' num2str(num_data_blocks)]); end;
+    clear buf raw
+    numvar = 1 + mod(k-1,noflds); %this is the numvar if storage is by variable
+    if strcmp(storvar,'var') & skip(numvar) == 1
+        fseek(fid,8192,0); %skip a block
+        continue
+    elseif lb == r8;
+        buf = fread(fid,lb,'double'); %read from real*8
+    else
+        %1365*5 = 6825
+        raw = char(fread(fid,6825,'5*uchar',1)); %read some real*5; skip the status byte
+        % The syntax of fread means you collect 6825 bytes of data in 5s, (ie 1365 elements) and 1365 skips of 1 byte.
+        % So there are 2 bytes left in the block of 8192.
+        fseek(fid,2,0); %skip the last 2 bytes
+        buf = punpack(raw(1:6825));
+    end
+    
+    %now store the data
+    
+    if strcmp(storvar,'var') %store this block in the correct variable
+        blocknum_thisvar = 1+floor((k-1)/noflds);
+        completed = lb*(blocknum_thisvar-1);
+        remaining = norecs-completed;
+        thistime = min(lb,remaining);
+        dcnums = completed+1:completed+thistime;
+        vn = varnames{numvar};
+        eval(['d.' vn '(dcnums) = buf(1:thistime);']);
+    else %store by data cycle; distribute data across variables
+        completed = numcycles_per_block*(k-1);
+        remaining = norecs-completed;
+        for k2 = 1:min(numcycles_per_block,remaining)
+            for k3 = loadlist
+                bufindex = noflds*(k2-1)+k3;
+                vn = varnames{k3};
+                eval(['d.' vn '(completed+k2) = buf(bufindex);']);
+            end %of k3 = 1:noflds
+        end %of k2 = 1:numcycles_per_block
+    end %of assigning data to variables
+
+end %of k = 1:num_data_blocks
+disp(['Finished reading data blocks']);
+
+fclose(fid);
+
+for k3 = loadlist  %sort out gridded data and absent values
+    clear grid vec absent bad
+    vn = varnames{k3};
+    eval(['vec = d.' vn ';']);
+    absent = h.absent(k3);
+    bad = find(vec == absent);
+    vec(bad) = nan;  %set absent data to nan;
+    if nrows > 1 & nplane < 2
+        vec = reshape(vec,nrows,ncols);
+    end
+    if nplane > 1
+        vec = reshape(vec,nrows,ncols,nplane);
+    end
+    eval(['d.' vn ' = vec;']);
+    %if nrows > 1 %assign data to 2-D matrix %better done with reshape command as above
+    %    for k1 = 1:nrows
+    %        for k2 = 1:ncols
+    %            grid(k1,k2) = vec(nrows*(k2-1)+k1);
+    %        end
+    %    end
+    %    eval(['d.' vn ' = grid;']); 
+    %else %copy one-D matrix with NaNs back to d.vn
+    %    eval(['d.' vn ' = vec;']);
+    %end
+end

pstar/progs/pload_old2.m

+function [d, h] = pload(filename,instring)
+
+% pload: load a pstar file; fixes variable names that would be illegal in matlab.
+% filename is a pstar filename
+% instring is an optional character array of variable numbers/names to load; use '/' or '' or ' ' for  all;
+% [d h] = pload(filename,instring);   %full form of command
+% [d h] = pload;              %will prompt for filename and varlist
+% [d h] = pload(filename);    %will prompt for variable list
+% [d h] = pload(filename,''); %will select all variables
+% [d h] = pload(filename,'1 3~5 a b');  %will select vars 1 3,4,5 and variables named 'a' and 'b'
+%
+% The final list of variable numbers is passed through 'unique' before loading.
+%
+% the variables are loaded into struct array d
+% the header info is loaded into struct array h
+%
+% BAK at SOC 31 March 2005
+
+% As presently written, the load is a bit slow for pstar data packed as
+% real*5; The time is mostly spent in punpack.m
+
+
+clear d h
+
+if nargin < 1
+    filename = input('Type name of pstar file to load :\n','s');
+elseif nargin > 2
+    disp(['You must give a maximum of two arguments']);
+    return
+end
+
+if exist(filename,'file') ~= 2
+    disp(['Filename ''' filename ''' not found']);
+    return
+end
+
+h = plisth(filename);
+norecs = h.norecs;
+noflds = h.noflds;
+nrows = max(h.nrows,1);
+%bak at NOCS, 14 Jul 2005: include possibility of NPLANE; I didn't know anyone actually used this,
+% but RTP says he sometimes does !
+nplane = max(h.nplane,1);
+
+if mod(norecs,nrows*nplane) ~= 0
+    error(['nrows * nplane ' num2str(nrows) ' * ' num2str(nplane) ' does not divide norecs ' num2str(norecs)])
+    return
+else
+    ncols = round(norecs/(nrows*nplane));
+end 
+
+if nargin < 2
+    disp([' ']);
+    disp(['Type list of pstar variables to load']);
+    disp(['List can be space or comma separated']);
+    instring = input('Use variable names or numbers or n1~n2. Type c/r or / for all); \n','s');
+end
+
+vlist = getvlist(instring,h);
+
+vlist = unique(vlist); %sort and remove duplicates for load
+
+skip = 1+0*[1:noflds];
+skip(vlist) = 0;
+loadlist = vlist;
+
+magvar=1344950870;
+magdat=1344950852;
+magvr8=1344951864;
+magdr8=1344947256;
+
+r5 = 1365;
+r8 = 1024;
+
+lb = 0;
+if h.magic == magvar | h.magic == magdat 
+    lb = r5;
+elseif h.magic == magvr8 | h.magic == magdr8
+    lb = r8;
+else
+    disp('Problem with magic number - quitting')
+    return
+end
+
+storvar = '   ';
+if h.magic == magvar | h.magic == magvr8
+    storvar = 'var';
+else
+    storvar = 'dcs';
+end
+
+varnames = '';
+disp(' ');
+printsummary = 0;  %Will set this to 1 if any variables are renamed
+for k = 1:noflds %sort out difficult variable names
+    if skip(k) == 1
+        varnames{k} = '';
+        continue
+    end
+    vn = h.fldnam{k};
+    vnold = vn;
+    renamed = 0;
+    %list of problem characters that may occur in pstar file names
+    %you can add to the list as new ones are discovered.
+    %Blanks are simply discarded.
+    if strcmp(vn,'        ') % if the pstar fldnam was all blanks, rename the variable.
+        vn = ['blank_' num2str(k) '_'];
+        renamed = 1;
+    end
+    
+    %remove trailing blanks;
+    while strcmp(vn(end),' ')
+        vn(end) = [];
+    end
+    
+    swap = {
+        ' ' '_space_'
+        '-' '_minus_'
+        '+' '_plus_'
+        '/' '_slash_'
+        '*' '_star_'
+        '.' '_dot_'
+        '#' '_hash_'
+        '$' '_dollar_'
+        '^' '_hat_'
+        '&' '_amp_'
+        '(' '_lparen_'
+        ')' '_rparen_'
+        '[' '_lbrac_'
+        ']' '_rbrac_'
+    } ;
+    nswaps = size(swap,1);
+    for k2 = 1:nswaps
+        s1 = char(swap(k2,1));
+        s2 = char(swap(k2,2));
+        sindex = findstr(vn,s1);
+        if length(sindex) > 0
+            vnew = '';
+            count = 1;
+            for kswap = 1:length(sindex)
+                vnew = [vnew vn(count:sindex(kswap)-1) s2];
+                count = sindex(kswap)+1;
+            end
+            vnew = [vnew vn(count:end)];
+            vn = vnew;
+            renamed = 1;
+        end
+    end
+    
+    %replace any other odd characters with underscore. We'll permit alphanumeric and underscore in matlab var names.
+    okchars = '_0123456789abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ';
+    for k2 = 1:length(vn)
+        if length(findstr(okchars,vn(k2))) == 0
+            vn(k2) = '_';
+            renamed = 1;
+        end
+    end
+
+    okfirstchar = 'abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ'; %must start with an alpha; if not, insert 'v_'
+    if length(findstr(okfirstchar,vn(1))) == 0
+        vn = ['v_' vn];
+        renamed = 1;
+    end
+    
+    if length(vn) > 26 %truncate. Allow for later addition of _nnn_ and still keep
+                       %total var name to shorter than 31, which seems to be needed in 
+                       % Matlab fieldnames in a structure.
+        vn = vn(1:26);
+        renamed = 1;
+    end
+    
+    
+    smatch = strmatch(vn,varnames,'exact');
+    if length(smatch) > 0
+        vn = [vn '_' num2str(k) '_'];
+        renamed = 1;
+    end
+    
+    if renamed == 1
+        printsummary = 1;
+        disp(['Variable number ' sprintf('%3d',k) ':''' vnold ''' renamed to ''' vn '''']);     
+    end
+    
+    vnprint = vn;
+    while length(vnprint) < 32
+        vnprint = [vnprint ' '];
+    end
+    varnames{k} = vnprint;
+end % of for k = 1:noflds %sort out problem variable names
+disp(' ');
+if printsummary == 1
+    disp('List of loaded variables will be');
+    disp('************************************************************************************************');
+    for k = loadlist
+        disp(['*' sprintf('%3d',k) '*' sprintf('%32s',varnames{k}) '*' h.fldunt{k} '* ' sprintf('%15.3f',h.alrlim(k)) ' * ' sprintf('%15.3f',h.uprlim(k)) ' * ' sprintf('%10.3f',h.absent(k)) ' *']);
+    end
+    disp('************************************************************************************************');
+    disp(' ')
+end
+
+if strcmp(storvar,'var') %store by variable
+    numblocks_per_var = ceil(norecs/lb);
+    num_data_blocks = noflds*numblocks_per_var;
+else %store by data cycle
+    numcycles_per_block = floor(lb/noflds);
+    num_data_blocks = ceil(norecs/numcycles_per_block);
+end
+
+disp(['There are ' num2str(num_data_blocks) ' data blocks to read'])
+
+fid = fopen(filename,'r','b');
+fseek(fid,8192,0); %skip header block
+
+for k = 1:num_data_blocks
+    if mod(k,100) == 0; disp(['Reading data block ' num2str(k) ' out of ' num2str(num_data_blocks)]); end;
+    clear buf raw
+    numvar = 1 + mod(k-1,noflds); %this is the numvar if storage is by variable
+    if strcmp(storvar,'var') & skip(numvar) == 1
+        fseek(fid,8192,0); %skip a block
+        continue
+    elseif lb == r8;
+        buf = fread(fid,lb,'double'); %read from real*8
+    else
+        %1365*5 = 6825
+        raw = char(fread(fid,6825,'5*uchar',1)); %read some real*5; skip the status byte
+        % The syntax of fread means you collect 6825 bytes of data in 5s, (ie 1365 elements) and 1365 skips of 1 byte.
+        % So there are 2 bytes left in the block of 8192.
+        fseek(fid,2,0); %skip the last 2 bytes
+        buf = punpack(raw(1:6825));
+    end
+    
+
+    %now store the data
+    if strcmp(storvar,'var') %store this block in the correct variable
+        blocknum_thisvar = 1+floor((k-1)/noflds);
+        completed = lb*(blocknum_thisvar-1);
+        remaining = norecs-completed;
+        thistime = min(lb,remaining);
+        dcnums = completed+1:completed+thistime;
+        vn = varnames{numvar};
+        if exist('d','var');
+            if isfield(d,m_remove_outside_spaces(vn))
+                eval(['lenvn = length(d.' vn ');']);
+                eval(['if lenvn < norecs; d.' vn '(end+1:norecs) = nan; end']);
+            end
+        end
+        eval(['d.' vn '(dcnums) = buf(1:thistime);']);
+
+    else %store by data cycle; distribute data across variables
+        completed = numcycles_per_block*(k-1);
+        remaining = norecs-completed;
+        for k2 = 1:min(numcycles_per_block,remaining)
+            for k3 = loadlist
+                bufindex = noflds*(k2-1)+k3;
+                vn = varnames{k3};
+                eval(['d.' vn '(completed+k2) = buf(bufindex);']);
+            end %of k3 = 1:noflds
+        end %of k2 = 1:numcycles_per_block
+    end %of assigning data to variables
+
+end %of k = 1:num_data_blocks
+disp(['Finished reading data blocks']);
+
+fclose(fid);
+
+for k3 = loadlist  %sort out gridded data and absent values
+    clear grid vec absent bad
+    vn = varnames{k3};
+    eval(['vec = d.' vn ';']);
+    absent = h.absent(k3);
+    bad = find(vec == absent);
+    vec(bad) = nan;  %set absent data to nan;
+    if nrows > 1 & nplane < 2
+        vec = reshape(vec,nrows,ncols);
+    end
+    if nplane > 1
+        vec = reshape(vec,nrows,ncols,nplane);
+    end
+    eval(['d.' vn ' = vec;']);
+    %if nrows > 1 %assign data to 2-D matrix %better done with reshape command as above
+    %    for k1 = 1:nrows
+    %        for k2 = 1:ncols
+    %            grid(k1,k2) = vec(nrows*(k2-1)+k1);
+    %        end
+    %    end
+    %    eval(['d.' vn ' = grid;']); 
+    %else %copy one-D matrix with NaNs back to d.vn
+    %    eval(['d.' vn ' = vec;']);
+    %end
+end

pstar/subs/getvlist.m

+function vlist = getvlist(list,h)
+
+% getvlist: compare character string 'list' with pstar fldnams in h.fldnam, and produce list of variable numbers
+% vlist = getvlist(list,h)
+%
+% list is comma or space delimited
+% list contains variable numbers or names, or a range defined by n~m,
+%    where n and m are variable numbers.
+% list has an optional trailing slash '/'
+% if list is empty, vlist is all variables.
+%
+% h is a pstar header
+%
+% BAK at SOC 31 March 2005
+
+vlist = [];
+for k = 1:h.noflds %remove blanks from var names
+    sname = h.fldnam{k};
+    i = strfind(sname,' ');
+    sname(i) = [];
+    fldnam{k} = sname;
+end
+
+% first swap comma to space
+i = strfind(list,',');
+while ~isempty(i)
+    list(i) = ' ';
+    i = strfind(list,',');
+end
+
+%remove trailing slash or space. Add trailing space first in case list
+%starts empty
+list = [list ' '];
+while strcmp(list(end),' ') | strcmp(list(end),'/')
+    list(end) = [];
+    if isempty(list); break; end
+end
+
+%remove all leading spaces
+list = [' ' list];
+while strcmp(list(1),' ')
+    list(1) = [];
+    if isempty(list); break; end
+end
+
+if isempty(list) % Empty variable list so return list of all variables
+    vlist = [1:h.noflds];
+    return
+end
+
+% remove multiple space
+i = strfind(list,'  ');
+while ~isempty(i)
+    list(i) = [];
+    i = strfind(list,'  ');
+end
+
+% The list should now be delimited by single spaces. Add one at each end
+list = [' ' list ' '];
+isp = strfind(list,' ');
+numv = length(isp)-1;
+
+badstr = 0;
+for k = 1:numv
+    vstr = list(isp(k)+1:isp(k+1)-1); %This is a string representing one of the elements of the input list
+    
+    i = findstr(vstr,'~'); %First look for range match
+    if isempty(i)  % No tilde
+    elseif length(i) > 1 % More than one tilde
+    elseif i == 1 | i == length(vstr) %tilde out of place
+    else %we have precisely one tilde, in a sensible place. Now look for valid variable numbers either side of it
+        vstr1 = vstr(1:i-1);
+        vstr2 = vstr(i+1:end);
+        t1 = testnum(vstr1,h.noflds);
+        t2 = testnum(vstr2,h.noflds);
+        if t1 == 1 & t2 == 1
+            vnum1 = str2num(vstr1);
+            vnum2 = str2num(vstr2);
+            vlist = [vlist vnum1:vnum2];
+            continue            
+        end
+    end
+        
+    % Now try for a simple variable number
+    t = testnum(vstr,h.noflds);
+    if t == 1
+        vnum = str2num(vstr);
+        vlist = [vlist vnum];
+        continue
+    end
+    
+    % Now try a string match
+    
+    vnum = [];
+    vnum = strmatch(vstr,fldnam,'exact');
+    if ~isempty(vnum) %variable string recognised; select the first if multiple matches.
+        vlist = [vlist vnum(1)];
+        continue
+    end
+    
+    % String not recognised
+    badstr = 1;
+    disp(['Variable string    ' vstr '    not recognised.']);
+end
+
+if badstr == 1
+    disp(' ');
+    error('One or more of the user-provided variable strings not recognised');
+end
+
+
+function t = testnum(str,noflds)
+
+% Test if the string represents a simple number in the range 1 to noflds
+% returns t = 1 if number is Ok. 0 otherwise.
+
+numerals = '0123456789';
+t = 1;
+for k = 1:length(str)
+    if isempty(findstr(numerals,str(k)))
+        t = 0;
+        return
+    end
+end
+
+num = str2num(str);
+if isempty(num)
+    t = 0;
+    return
+end
+
+if round(num) ~= num
+    t = 0;
+    return
+end
+
+if num < 1 | num > noflds
+    t = 0;
+    return
+end

pstar/subs/punpack.m

+punpack_old.m
\ No newline at end of file

pstar/subs/punpack2.m

+function out = punpack(in)
+
+% punpack: unpack real*5 data in pstar files
+% mfile to unpack real*5 data for loading pstar files.
+% Use:
+%  out = punpack(in)
+%
+% The array in contains a series of 5-byte
+% representations of real numbers, read in as characters char*5.
+% The pstar status byte was skipped at read time, to make this subroutine
+%  more efficient, and enable use of simple reshape.
+%
+% These char*5 are converted to 5 times hex*2, and then concatenated and converted to real*8.
+%
+%  Note that hex2num pads from real*5 to real*8 with zeros, as required for pstar real*5 files.
+%
+% BAK at SOC 31 March 2005
+
+%much faster version by BAK 23 Jun 2007
+%exploit knowledge of ieee real numbers, instead of matlab call to HEX2NUM,
+% as in previous version of punpack.
+
+in = uint8(in);  %force to uint8
+l = length(in)/5;
+in = reshape(in,5,l);
+in = in';
+a = double(in);  %perform arithmetic with double, so that sign is carried
+
+% a real*8 consists of 1 sign bit; 11 exponent bits; 52 fraction bits
+% a pstar real*5 is the above with the last 3 bytes discarded
+
+
+
+% first, split byte 2 of each real*5 into 4 most and least significant bits
+
+b2 = in(:,2);
+b2_l = double((b2-8)/16); %innermost divide is with uint8; 
+% this truncates the 4 right-hand bits in the byte
+% leaving the 4 left hand bits in the range 0000 to 1111
+b2_r = a(:,2) - 16*b2_l; %now recover the 4 right hand bits
+
+
+a1 = in(:,1);
+neg = find(a1 > 127);        %find where sign bit set
+a(neg,1) = a(neg,1)-128;     %adjust for sign bit in first byte
+
+% exponent is now first 12 bits, allowing us to include adjusted sign bit
+
+e = a(:,1)*16 + b2_l - 1023;  %exponent is stored offset by 1023
+
+f = 1 + (((a(:,5)/256 + a(:,4))/256 + a(:,3))/256 + b2_r)/16; %build fraction
+
+out = pow2(f,e);
+out(neg) = -out(neg);
+return
\ No newline at end of file

pstar/subs/punpack3.m

+function out = punpack(in)
+% disp('p3')
+
+%much faster version by BAK 23 Jun 2007
+%exploit knowledge of ieee real numbers, instead of matlab call to HEX2NUM,
+% as in previous version of punpack.
+
+in = uint8(in);  %force to uint8
+l = length(in)/5;
+in = reshape(in,5,l);
+in = in';
+a = double(in);  %perform arithmetic with double, so that sign is carried
+
+% a real * 8 consists of 1 sign bit; 11 exponent bits; 52 fraction bits
+% a pstar real*5 is the above with the last 3 bytes discarded
+
+% split array of second bytes in real*5 into 4 most and least significant
+% bits
+
+a2_l = floor(a(:,2)/16); % get the most significant part, equivalent to the 4 left hand bits
+a2_r = a(:,2) - 16*a2_l; % now recover the 4 right hand bits which are part of the fraction
+
+a1 = in(:,1);
+neg = find(a1 > 127);        %find where sign bit set
+a(neg,1) = a(neg,1)-128;     %adjust for sign bit in first byte
+
+% exponent is now first 12 bits, allowing us to include adjusted sign bit
+
+e = a(:,1)*16 + a2_l - 1023;  %exponent is stored offset by 1023
+
+f = 1 + (((a(:,5)/256 + a(:,4))/256 + a(:,3))/256 + a2_r)/16; %build fraction
+
+out = pow2(f,e);
+out(neg) = -out(neg);
+return
\ No newline at end of file

pstar/subs/punpack4.m

+function out = punpack(in)
+% disp('p3')
+
+fid = fopen('/users/bak/wkk','r+');
+
+
+%much faster version by BAK 23 Jun 2007
+%exploit knowledge of ieee real numbers, instead of matlab call to HEX2NUM,
+% as in previous version of punpack.
+
+in = uint8(in);  %force to uint8
+l = length(in)/5;
+
+in = in(:)';
+
+z = uint8(0);
+kk = 0;
+o = [];
+for k = 1:l
+    oo = [in(kk+1:kk+5) z z z];
+    kk = kk+5;
+    o = [o oo];
+end
+fwrite(fid,o,'uint8');
+fseek(fid,0,-1);
+out = fread(fid,l,'double');
+% if l > 1 ; keyboard; end
+fclose(fid);
+
+return

pstar/subs/punpack5.m

+function out = punpack(in)
+
+% punpack: unpack real*5 data in pstar files
+% mfile to unpack real*5 data for loading pstar files.
+% Use:
+%  out = punpack(in)
+%
+% The array in contains a series of 5-byte
+% representations of real numbers, read in as characters char*5.
+% The pstar status byte was skipped at read time, to make this subroutine
+%  more efficient, and enable use of simple reshape.
+%
+% These char*5 are converted to 5 times hex*2, and then concatenated and converted to real*8.
+%
+%  Note that hex2num pads from real*5 to real*8 with zeros, as required for pstar real*5 files.
+%
+% BAK at SOC 31 March 2005
+
+%much faster version by BAK 23 Jun 2007
+%exploit knowledge of ieee real numbers, instead of matlab call to HEX2NUM,
+% as in previous version of punpack.
+
+% if ~strcmp(class(in),'double'); in = uint8(in);  %force to uint8
+l = length(in)/5;
+in = reshape(in,5,l);
+in = in';
+if ~strcmp(class(in),'double');in = double(in); end  %force to uint8
+% a = double(in);  %perform arithmetic with double, so that sign is carried
+
+% a real*8 consists of 1 sign bit; 11 exponent bits; 52 fraction bits
+% a pstar real*5 is the above with the last 3 bytes discarded
+
+
+
+% first, split byte 2 of each real*5 into 4 most and least significant bits
+
+% % % b2 = in(:,2);
+% % % b2_l = double((b2-8)/16); %innermost divide is with uint8; 
+% % % % this truncates the 4 right-hand bits in the byte
+% % % % leaving the 4 left hand bits in the range 0000 to 1111
+% % % b2_r = a(:,2) - 16*b2_l; %now recover the 4 right hand bits
+
+b2_l = floor(in(:,2)/16);
+b2_r = in(:,2)-16*b2_l;
+
+a1 = in(:,1);
+neg = find(a1 > 127);        %find where sign bit set
+in(neg,1) = in(neg,1)-128;     %adjust for sign bit in first byte
+
+% exponent is now first 12 bits, allowing us to include adjusted sign bit
+
+e = in(:,1)*16 + b2_l - 1023;  %exponent is stored offset by 1023
+
+f = 1 + (((in(:,5)/256 + in(:,4))/256 + in(:,3))/256 + b2_r)/16; %build fraction
+
+out = pow2(f,e);
+out(neg) = -out(neg);
+return
\ No newline at end of file

pstar/subs/punpack_old.m

+function out = punpack(in)
+
+% punpack: unpack real*5 data in pstar files
+% mfile to unpack real*5 data for loading pstar files.
+% Use:
+%  out = punpack(in)
+%
+% The array in contains a series of 5-byte
+% representations of real numbers, read in as characters char*5.
+% The pstar status byte was skipped at read time, to make this subroutine
+%  more efficient, and enable use of simple reshape.
+%
+% These char*5 are converted to 5 times hex*2, and then concatenated and converted to real*8.
+%
+%  Note that hex2num pads from real*5 to real*8 with zeros, as required for pstar real*5 files.
+%
+% BAK at SOC 31 March 2005
+
+n = length(in); %array 'in' should be an exact multiple of 5 bytes.
+
+call = sprintf('%02x',in);
+
+c = reshape(call,10,n/5)';
+
+out = hex2num(c)';
+
+return

pstar/subs/punpackn.m

+function out = punpack(in)
+
+%much faster version by BAK 23 Jun 2007
+%exploit knowledge of ieee real numbers, instead of matlab call to HEX2NUM,
+% as in previous version of punpack.
+
+in = uint8(in);  %force to uint8
+l = length(in)/5;
+in = reshape(in,5,l);
+in = in';
+a = double(in);  %perform arithmetic with double, so that sign is carried
+
+% a real * 8 consists of 1 sign bit; 11 exponent bits; 52 fraction bits
+% a pstar real*5 is the above with the last 3 bytes discarded
+
+% split array of second bytes in real*5 into 4 most and least significant
+% bits
+
+a2_l = floor(a(:,2)/16); % get the most significant part, equivalent to the 4 left hand bits
+a2_r = a(:,2) - 16*a2_l; % now recover the 4 right hand bits which are part of the fraction
+
+a1 = in(:,1);
+neg = find(a1 > 127);        %find where sign bit set
+a(neg,1) = a(neg,1)-128;     %adjust for sign bit in first byte
+
+% exponent is now first 12 bits, allowing us to include adjusted sign bit
+
+e = a(:,1)*16 + a2_l - 1023;  %exponent is stored offset by 1023
+
+f = 1 + (((a(:,5)/256 + a(:,4))/256 + a(:,3))/256 + a2_r)/16; %build fraction
+
+out = pow2(f,e);
+out(neg) = -out(neg);
+return
\ No newline at end of file