accome_yp

 This function reads the equilibrium parameters from file created by ACCOME and saved in the formatted text file "filename"

 by Y.PEYSSON (CEA-DRFC, yves.peysson@cea.fr)

0001 function [x,y,xypsin,psia,rpsin,xyBPHI,xyBR,xyBZ,Rp,Zp] = accome_yp(filename,display_mode,p_opt,dispname,savename,Ns);
0002 %
0003 % This function reads the equilibrium parameters from file created by ACCOME and saved in the formatted text file "filename"
0004 %
0005 % by Y.PEYSSON (CEA-DRFC, yves.peysson@cea.fr)
0006 %
0007 if nargin < 6
0008     Ns = 1001;%magnetic axis recalculation supergrid factor
0009 end
0010 if nargin < 5
0011     savename = filename;
0012 end
0013 if nargin < 4
0014     dispname = filename;
0015 end
0016 %
0017 fid = fopen(filename);
0018 %
0019 fscanf(fid,'%s',5);
0020 Rp = fscanf(fid,'%f',1);%major radius of magnetic axis
0021 Zp = fscanf(fid,'%f',1);%vertical height of magnetic axis
0022 %
0023 fscanf(fid,'%s',6);
0024 nxminmax = fscanf(fid,'%d',2);nx = nxminmax(2) - nxminmax(1) + 1;
0025 nyminmax = fscanf(fid,'%d',2);ny = nyminmax(2) - nyminmax(1) + 1;
0026 %
0027 fscanf(fid,'%s',2);
0028 R = fscanf(fid,'%f',nx);
0029 x = R - Rp;
0030 fscanf(fid,'%s',2);
0031 Z = fscanf(fid,'%f',ny);
0032 y = Z - Zp;
0033 %
0034 fscanf(fid,'%s',1);
0035 xyBPHI = fscanf(fid,'%f',nx*ny);%Toroidal magnetic field (T)
0036 xyBPHI = reshape(xyBPHI,nx,ny);
0037 %
0038 fscanf(fid,'%s',3);
0039 xypsi = fscanf(fid,'%f',nx*ny);%Poloidal magnetic flux function
0040 xypsi = reshape(xypsi,nx,ny);
0041 %
0042 fscanf(fid,'%s',1);
0043 xyBR = fscanf(fid,'%f',nx*ny);%Poroidal magnetic field BR (T)
0044 if ~isempty(xyBR),
0045     xyBR = reshape(xyBR,nx,ny);
0046 end
0047 %
0048 fscanf(fid,'%s',1);
0049 xyBZ = fscanf(fid,'%f',nx*ny);%Poroidal magnetic field BZ (T)
0050 if ~isempty(xyBZ),
0051     xyBZ = reshape(xyBZ,nx,ny);
0052 end
0053 %
0054 fclose(fid);
0055 %
0056 %Automatic search of last closed magnetic flux surface and the magnetic
0057 %flux values
0058 %
0059 psic = contourc(x,y,xypsi',Ns);
0060 %
0061 %figure(100),contour(x,y,xypsi',Ns);axis('equal');
0062 %
0063 psi_offset = 1;
0064 for ic = 1:Ns
0065     psi_level(ic) = psic(1,psi_offset);
0066     psi_index(ic) = psic(2,psi_offset);
0067 %
0068 %    psi_v = psic(:,psi_offset+1:psi_offset + psi_index(ic));
0069 %    if (psi_v(1,1) == psi_v(1,end)) & (psi_v(2,1) == psi_v(2,end)), psi_status(ic) = 1;else psi_status(ic) = 0;end
0070 %
0071     psi_offset = psi_index(ic) + psi_offset + 1;
0072 end
0073 %
0074 psip = psi_level(1);%poloidal flux function value at the magnetic axis
0075 %psia = psi_level(find(psi_status==0,1)-1);%poloidal flux function value at the last closed flux surface
0076 psia = psi_level(find(psi_level>=0,1)-1);%poloidal flux function value at the last closed flux surface
0077 %
0078 xypsi = xypsi - psip;%psi = 0 on the magnetic axis (arbitrary choice)
0079 psia = (psia - psip);%if > 0, IPHI is positive, PHI clockwise from top otherwise IPHI is negative
0080 %
0081 xypsin = xypsi/psia;
0082 rpsin = linspace(0,1,nx);
0083 %
0084 save(['test_EQUIL_',savename,'.mat'])
0085