file naming for WDM simulation

projects/WDM/WDM_model.m

@@ -1,25 +1,96 @@
 %%% Run parameters
 % TX
-% --- FIRST LINE: evaluate settings located beside this script ---
-run(fullfile(fileparts(mfilename('fullpath')),'WDM_settings.m'));
 
-num_realiz      = 50;
+function WDM_model(options)
-s.wavelengthplan  = calcWavelengthPlan(16,400e9,1310);
+
+arguments
+    options.num_channels = 16;
+    options.channel_spacing = 400e9;
+    options.fiber_length_km = 0;
+    options.rand_key = 1;
+    options.num_realiz = 1;
+    options.fwm_mitigation_technique = "co";
+end
+
+%% 
+
+% Add the imdd_simulation framework to the path
+if ispc
+    addpath(genpath('C:\Users\Silas\Documents\MATLAB\imdd_simulation'));
+else
+    % Linux path on the cluster
+    addpath(genpath('/work_beegfs/sutef391/imdd_simulation'));
+end
+
+% Quiet the ambiguous CET warning (best is to set TZ in sbatch; see below)
+warning('off','MATLAB:datetime:AmbiguousTimeZone');
+
+% How many workers?
+cpus = str2double(getenv('SLURM_CPUS_PER_TASK'));
+if ~isfinite(cpus) || cpus < 1, cpus = max(1, feature('numcores')); end
+
+% Use a per-job, node-local JobStorageLocation to avoid stale locks on $HOME
+% Prefer $TMPDIR if your cluster provides it, else tempdir().
+tmpbase = getenv('TMPDIR');
+if isempty(tmpbase), tmpbase = tempdir; end
+jsl = fullfile(tmpbase, sprintf('matlab_jobstorage_%s_%s', ...
+    getenv('USER'), getenv('SLURM_JOB_ID')));
+if ~exist(jsl,'dir'); mkdir(jsl); end
+
+% Configure the local cluster explicitly and start the pool
+c = parcluster('local');
+c.NumWorkers = cpus;
+c.JobStorageLocation = jsl;
+
+p = gcp('nocreate');
+if isempty(p) || p.NumWorkers ~= cpus
+    if ~isempty(p), delete(p); end
+    p = parpool(c, cpus);           % avoids the “queued” state
+end
+fprintf('parpool up with %d workers; JobStorage=%s\n', p.NumWorkers, c.JobStorageLocation);
+
+% result filename (timestamp + optional job id)
+t      = datetime('now','TimeZone','local','Format','yyyyMMdd_HHmmss');
+jobid  = getenv('SLURM_JOB_ID'); if isempty(jobid), jobid = 'nojid'; end
+host   = getenv('HOSTNAME'); if isempty(host), host = 'localhost'; end
+
+% Output directory depends on platform
+% create/ use folders foroptions.fiber_length_km, options.num_channels, options.channel_spacing.*1e-9, options.fwm_mitigation_technique
+foldname = sprintf('%dkm_%dch_%dghz_%s', options.fiber_length_km, options.num_channels, options.channel_spacing.*1e-9, options.fwm_mitigation_technique);
+if ispc
+    output_root = fullfile('C:\Users\Silas\Documents\MATLAB\Datensätze\FWM_2025\',foldname,'\');
+else
+    output_root =  fullfile('/work_beegfs/sutef391/results_WDM',foldname,'\');
+end
+if ~exist(output_root,'dir'), mkdir(output_root); end
+
+% Build filename
+t     = datetime('now','TimeZone','local','Format','yyyyMMdd_HHmmss');
+jobid = getenv('SLURM_JOB_ID'); if isempty(jobid), jobid = 'nojid'; end
+host  = getenv('HOSTNAME');     if isempty(host),  host  = 'localhost'; end
+
+fname = sprintf('WDM_%s_%s_%s_%dkm_%dch_%dghz_%s.mat', char(t), host, jobid, options.fiber_length_km, options.num_channels, options.channel_spacing.*1e-9, options.fwm_mitigation_technique);
+
+
+%%
+s.num_realiz      = options.num_realiz;
+% s.wavelengthplan  = calcWavelengthPlan(16,400e9,1310);
+s.wavelengthplan  = calcWavelengthPlan(options.num_channels,options.channel_spacing,1310);
 % wavelengthplan  = [1295,1305,1315,1325];
-link_length     = 10;
+link_length       = options.fiber_length_km;
-pmd             = 0.1;
+s.pmd             = 0.1;
-gamma           = 0.0023;
+s.gamma           = 0.0023;
 
-M               = 4;
+s.M               = 4;
-m = floor(log2(M)*10)/10;
+m = floor(log2(s.M)*10)/10;
 fsym            = 112e9;
 fdac            = 2*fsym;
-fadc            = 2*fsym;
+fadc            = 120000000000;
-s.random_key      = 100;
+s.random_key      = options.rand_key;
 
 % Laser / s.Modulator
 vbias_rel       = 0.5;
-u_pi            = 3.2;
+u_pi            = 4.6;
 vbias           = -vbias_rel*u_pi;
 laser_linewidth = 0e6;
 
@@ -52,7 +123,7 @@ f_plan = physconst('lightspeed')./(s.wavelengthplan.*1e-9);
 margin = 25e12; % some THz left and right
 f_span = (max(f_plan)+margin)-(min(f_plan)-margin);
 f_nyq = f_span/2;
-kover = 8;
+kover = 4;
 upsample_required = f_nyq./(fdac*kover/2);
 upsample_pow = 2^nextpow2(upsample_required);
 upsample_ceil = ceil(upsample_required);
@@ -64,13 +135,29 @@ signal_cell = {};
 Symbols = {};
 Tx_bits = {};
 
-s.rop = -6:0.75:-0.75;
+s.rop = -12:0.75:-0.75;
 
 output_ffe = cell(length(s.wavelengthplan),length(s.rop),s.num_realiz);
 output_vnle = cell(length(s.wavelengthplan),length(s.rop),s.num_realiz);
 output_mlse = cell(length(s.wavelengthplan),length(s.rop),s.num_realiz);
 output_dbt = cell(length(s.wavelengthplan),length(s.rop),s.num_realiz);
 
+s.p = options.fwm_mitigation_technique;
+switch s.p
+    case "co"
+        pol_rot = 100.*ones(1,length(s.wavelengthplan));
+        d_local = 0;
+    case "pair"
+        pol_rot = repmat([100,100,0,0],1,length(s.wavelengthplan)/4);
+        d_local = 0;
+    case "alt"
+        pol_rot = repmat([100,0,100,0],1,length(s.wavelengthplan)/4);
+        d_local = 0;
+    case "seg"
+        pol_rot = 100.*ones(1,length(s.wavelengthplan));
+        d_local = 3;
+end
+
 for realiz = 1:s.num_realiz
 
 
@@ -86,24 +173,25 @@ for realiz = 1:s.num_realiz
             "mrds_code",0,"mrds_blocklength",512,"duobinary_mode",duob_mode).process();
 
         % Digi_sig.spectrum("fignum",101,"displayname",'bla','normalizeTo0dB',0,'lambda0_nm',1310,'useWavelengthAxis',0);
-        Lp_awg = Filter('filtdegree',3,"f_cutoff",100e9,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true);
+        Lp_awg = Filter('filtdegree',3,"f_cutoff",56e9,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true);
-        El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",1,"kover",kover,"bit_resolution",12,"upsampling_method","samplehold","precomp_sinc_rolloff",0,"H_lpf",Lp_awg,"dac_max",0.6,"dac_min",-0.6).process(Digi_sig);
+        El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",1,"kover",kover,"bit_resolution",6,"upsampling_method","samplehold","precomp_sinc_rolloff",0,"H_lpf",Lp_awg,"dac_max",0.6,"dac_min",-0.6).process(Digi_sig);
         % El_sig = s.M8199B("kover",kover).process(Digi_sig);
         % El_sig.spectrum("fignum",101,"displayname",'bla','normalizeTo0dB',0,'lambda0_nm',1310,'useWavelengthAxis',0);
 
         %%%%% Electrical Driver Amplifier %%%%%%
         El_sig = El_sig.normalize("mode","oneone");
+        El_sig = El_sig .* u_pi .* 0.5;
         % El_sig = El_sig.setPower(1,"dBm");
         % figure;histogram(El_sig.signal);
 
         %%%%% s.MODULATE E/O CONVERSION %%%%%
         Eml_out = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",s.wavelengthplan(l),"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",s.random_key+l+realiz).process(El_sig);
 
-        signal_cell{l} = Polarization_Controller("mode","rot_power","desired_power",100).process(Eml_out);
+        signal_cell{l} = Polarization_Controller("mode","rot_power","desired_power",pol_rot(l)).process(Eml_out);
     end
 
     Opt_sig_wdm = Optical_Multiplex("fs_in",fdac*kover,"fs_out",upsample_pow*fdac*kover,...
-                                "lambda_center",1310,"random_key",0,"filtype",1,"B",200e9).process(signal_cell);
+        "lambda_center",1310,"random_key",0,"filtype",1,"B",120e9).process(signal_cell);
 
     Opt_sig_wdm = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",3+10*log10(N)).process(Opt_sig_wdm);
 
@@ -117,37 +205,35 @@ for realiz = 1:s.num_realiz
     segment_length = 1; % km
     nSegments = link_length/segment_length;
     zdw = 1310;
-    D_local = 0; %if ~=0, simulation uses "segmented fiber with d+,d-)
+    d_local = d_local; %if ~=0, simulation uses "segmented fiber with d+,d-)
     randomize_D = true;
-    Dvec = getDispersionVector(nSegments, D_local, zdw, randomize_D, s.random_key+realiz);
+    Dvec = getDispersionVector(nSegments, d_local, zdw, randomize_D, s.random_key+realiz);
     for seg = 1:nSegments
 
-        Opt_sig_wdm_fib = DP_Fiber("L",segment_length,"D",Dvec(s),"Dpmd",pmd,"Ds",0.07,...
+        Opt_sig_wdm_fib = DP_Fiber("L",segment_length,"D",Dvec(seg),"Dpmd",s.pmd,"Ds",0.07,...
             "beat_len",10,"corr_len",100,"dz",1,"manakov",0,...
             "gamma",s.gamma,"lambda",zdw,"n_waveplates",10,"SS_dphimax",0.01,...
             "SS_dzmax",50,"SS_dzmin",10,"X_alpha",0.3,"X_beta",0,"rng",1).process(Opt_sig_wdm_fib);
 
     end
 
-    Opt_sig_wdm_fib.spectrum("fignum",realiz,"displayname",'bla','lambda0_nm',1310,'useWavelengthAxis',0);
+    % Opt_sig_wdm_fib.spectrum("fignum",realiz,"displayname",'bla','lambda0_nm',1310,'useWavelengthAxis',0);
 
     % Opt_sig_wdm_fib.move_it_spectrum("fignum",100212,"displayname",'bla');
 
     % Opt_sig    = Fiber("fsimu",Opt_sig.fs,"fiber_length",s.link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"s.gamma",0,"Dslope",0.07).process(Opt_sig)
+    Opt_sig_wdm_demux = Optical_Demultiplex("attenuation",0,"B",200e9,"filtype",1,"fs_out",fdac*kover,"fs_in",fdac*kover*upsample_pow,"lambda_center",1310).process(Opt_sig_wdm_fib);
 
     for ri = 1:length(s.rop)
 
+        parfor l = 1:N
+
             %%%%%% ROP %%%%%%
-        Opt_sig_wdm_rx = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",s.rop(ri)+10*log10(N)).process(Opt_sig_wdm_fib);
+            Opt_sig_wdm_rx = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",s.rop(ri)).process(Opt_sig_wdm_demux{l}); % rop+10*log10(N)
-        
-        Opt_sig_wdm_demux = Optical_Demultiplex("attenuation",0,"B",200e9,"filtype",1,"fs_out",Opt_sig_wdm_rx.fs/upsample_pow,"fs_in",Opt_sig_wdm_rx.fs,"lambda_center",1310).process(Opt_sig_wdm_rx);
-        
-        PD_cell = {};
-        for l = 1:N
 
             %%%%%% PD Square Law %%%%%%
-            assert(fdac*kover==Opt_sig_wdm_demux{l}.fs,'Sampling Frequencies do not match! Check previous steps');
+            assert(fdac*kover==Opt_sig_wdm_rx.fs,'Sampling Frequencies do not match! Check previous steps');
-            PD_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20,"nep",1.8e-11,"randomkey",s.random_key+l+realiz).process(Opt_sig_wdm_demux{l});
+            PD_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20,"nep",1.8e-11,"randomkey",s.random_key+l+realiz).process(Opt_sig_wdm_rx);
 
             % PD_sig.spectrum("fignum",222,"displayname",'bla','normalizeTo0dB',1);
 
@@ -156,13 +242,13 @@ for realiz = 1:s.num_realiz
             PD_sig = Filter('filtdegree',4,"f_cutoff",rx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(PD_sig);
 
             % %%%%%% Low-pass Scope %%%%%%
-            Lp_scpe = Filter('filtdegree',4,"f_cutoff",110e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true);
+            Lp_scpe = Filter('filtdegree',4,"f_cutoff",80e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true);
 
             %%%%%% Scope %%%%%%
             Scpe_sig = Scope("fsimu",fdac*kover,"fadc",fadc,...
                 "delay",0,"fixed_delay",0,"filtertype",filtertypes.butterworth,...
                 "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,...
-                "adcresolution",8,"quantbuffer",0.1,'block_dc',1,'lpf_active',1,'H_lpf',Lp_scpe).process(PD_sig);
+                "adcresolution",8,"quantbuffer",0.1,'block_dc',1,'lpf_active',0,'H_lpf',Lp_scpe).process(PD_sig);
 
             Scpe_sig_2sps = Scpe_sig.resample("fs_out",2*fsym);
             % Scpe_sig.spectrum("fignum",222,"displayname",'bla','normalizeTo0dB',1);
@@ -244,36 +330,7 @@ for realiz = 1:s.num_realiz
 
 
 end
-
-function dispersion_vector = getDispersionVector(N, D, ref_zdw, randomize_ZDW, randomkey)
-% s.MATLAB version of the Python generator shown above.
-% Returns an N×1 vector (ps/(nm·km)).
-%
-% D is the nominal dispersion magnitude. For D>0 the link is segmented with
-% alternating sign (+D, -D, +D, …). For D==0 it is flat (0) except for
-% ZDW randomization. The ZDW detuning is ~N(0, 2 nm) around 1310 nm and is
-% converted to dispersion via 0.09 ps/(nm·km) per nm.
-
-    % constants (matching the Python code)
-    meanLambda_nm = 1310;        % center wavelength
-    sigma_nm      = 2;           % ZDW sigma
-    Dslope        = 0.07;        % ps/(nm·km) per nm detuning
-
-    % random ZDW-induced dispersion offset
-    if randomize_ZDW
-        rng(randomkey, 'twister');
-        rand_zdws_nm = meanLambda_nm + sigma_nm .* randn(N,1);
-        rand_D = (rand_zdws_nm - ref_zdw) .* Dslope;   % ps/(nm·km)
-    else
-        rand_D = zeros(N,1);
 end
 
-    % nominal segmented pattern (match Python intent; keep length N)
-    if D > 0
-        base = (-1) .^ ((0:N-1).');   % +1,-1,+1,-1,...
-    else % D == 0 (or anything else)
-        base = ones(N,1);
-    end
 
-    dispersion_vector = base .* D + rand_D;  % ps/(nm·km)
+
-end

projects/WDM/WDM_model_old.m

@@ -3,23 +3,23 @@
 % --- FIRST LINE: evaluate settings located beside this script ---
 run(fullfile(fileparts(mfilename('fullpath')),'WDM_settings.m'));
 
-s.num_realiz      = 2;
+num_realiz      = 50;
 s.wavelengthplan  = calcWavelengthPlan(16,400e9,1310);
 % wavelengthplan  = [1295,1305,1315,1325];
-link_length     = 0;
+link_length     = 10;
 pmd             = 0.1;
 gamma           = 0.0023;
 
-s.M               = 4;
+M               = 4;
-m = floor(log2(s.M)*10)/10;
+m = floor(log2(M)*10)/10;
 fsym            = 112e9;
 fdac            = 2*fsym;
-fadc            = 120000000000;
+fadc            = 2*fsym;
 s.random_key      = 100;
 
 % Laser / s.Modulator
 vbias_rel       = 0.5;
-u_pi            = 4.6;
+u_pi            = 3.2;
 vbias           = -vbias_rel*u_pi;
 laser_linewidth = 0e6;
 
@@ -86,14 +86,13 @@ for realiz = 1:s.num_realiz
             "mrds_code",0,"mrds_blocklength",512,"duobinary_mode",duob_mode).process();
     
         % Digi_sig.spectrum("fignum",101,"displayname",'bla','normalizeTo0dB',0,'lambda0_nm',1310,'useWavelengthAxis',0);
-        Lp_awg = Filter('filtdegree',3,"f_cutoff",56e9,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true);
+        Lp_awg = Filter('filtdegree',3,"f_cutoff",100e9,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true);
-        El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",1,"kover",kover,"bit_resolution",6,"upsampling_method","samplehold","precomp_sinc_rolloff",0,"H_lpf",Lp_awg,"dac_max",0.6,"dac_min",-0.6).process(Digi_sig);
+        El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",1,"kover",kover,"bit_resolution",12,"upsampling_method","samplehold","precomp_sinc_rolloff",0,"H_lpf",Lp_awg,"dac_max",0.6,"dac_min",-0.6).process(Digi_sig);
         % El_sig = s.M8199B("kover",kover).process(Digi_sig);
         % El_sig.spectrum("fignum",101,"displayname",'bla','normalizeTo0dB',0,'lambda0_nm',1310,'useWavelengthAxis',0);
     
         %%%%% Electrical Driver Amplifier %%%%%%
         El_sig = El_sig.normalize("mode","oneone");
-        El_sig = El_sig .* u_pi .* 0.5;
         % El_sig = El_sig.setPower(1,"dBm");
         % figure;histogram(El_sig.signal);
     
@@ -104,7 +103,7 @@ for realiz = 1:s.num_realiz
     end
     
     Opt_sig_wdm = Optical_Multiplex("fs_in",fdac*kover,"fs_out",upsample_pow*fdac*kover,...
-                                "lambda_center",1310,"random_key",0,"filtype",1,"B",120e9).process(signal_cell);
+                                "lambda_center",1310,"random_key",0,"filtype",1,"B",200e9).process(signal_cell);
     
     Opt_sig_wdm = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",3+10*log10(N)).process(Opt_sig_wdm);
     
@@ -135,18 +134,20 @@ for realiz = 1:s.num_realiz
     % Opt_sig_wdm_fib.move_it_spectrum("fignum",100212,"displayname",'bla');
     
     % Opt_sig    = Fiber("fsimu",Opt_sig.fs,"fiber_length",s.link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"s.gamma",0,"Dslope",0.07).process(Opt_sig)
-    Opt_sig_wdm_demux = Optical_Demultiplex("attenuation",0,"B",200e9,"filtype",1,"fs_out",Opt_sig_wdm_rx.fs/upsample_pow,"fs_in",Opt_sig_wdm_rx.fs,"lambda_center",1310).process(Opt_sig_wdm_fib);
     
     for ri = 1:length(s.rop)
     
-        parfor l = 1:N
-
         %%%%%% ROP %%%%%%
-            Opt_sig_wdm_rx = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",s.rop(ri)).process(Opt_sig_wdm_demux{l}); % rop+10*log10(N)
+        Opt_sig_wdm_rx = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",s.rop(ri)+10*log10(N)).process(Opt_sig_wdm_fib);
+        
+        Opt_sig_wdm_demux = Optical_Demultiplex("attenuation",0,"B",200e9,"filtype",1,"fs_out",Opt_sig_wdm_rx.fs/upsample_pow,"fs_in",Opt_sig_wdm_rx.fs,"lambda_center",1310).process(Opt_sig_wdm_rx);
+        
+        PD_cell = {};
+        for l = 1:N
         
             %%%%%% PD Square Law %%%%%%
-            assert(fdac*kover==Opt_sig_wdm_rx.fs,'Sampling Frequencies do not match! Check previous steps');
+            assert(fdac*kover==Opt_sig_wdm_demux{l}.fs,'Sampling Frequencies do not match! Check previous steps');
-            PD_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20,"nep",1.8e-11,"randomkey",s.random_key+l+realiz).process(Opt_sig_wdm_rx);
+            PD_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20,"nep",1.8e-11,"randomkey",s.random_key+l+realiz).process(Opt_sig_wdm_demux{l});
         
             % PD_sig.spectrum("fignum",222,"displayname",'bla','normalizeTo0dB',1);
         
@@ -155,13 +156,13 @@ for realiz = 1:s.num_realiz
             PD_sig = Filter('filtdegree',4,"f_cutoff",rx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(PD_sig);
         
             % %%%%%% Low-pass Scope %%%%%%
-            Lp_scpe = Filter('filtdegree',4,"f_cutoff",80e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true);
+            Lp_scpe = Filter('filtdegree',4,"f_cutoff",110e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true);
         
             %%%%%% Scope %%%%%%
             Scpe_sig = Scope("fsimu",fdac*kover,"fadc",fadc,...
                 "delay",0,"fixed_delay",0,"filtertype",filtertypes.butterworth,...
                 "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,...
-                "adcresolution",8,"quantbuffer",0.1,'block_dc',1,'lpf_active',0,'H_lpf',Lp_scpe).process(PD_sig);
+                "adcresolution",8,"quantbuffer",0.1,'block_dc',1,'lpf_active',1,'H_lpf',Lp_scpe).process(PD_sig);
         
             Scpe_sig_2sps = Scpe_sig.resample("fs_out",2*fsym);
             % Scpe_sig.spectrum("fignum",222,"displayname",'bla','normalizeTo0dB',1);

projects/WDM/WDM_settings.m

@@ -2,60 +2,3 @@
 
 
 
-% Add the imdd_simulation framework to the path
-if ispc
-    addpath(genpath('C:\Users\Silas\Documents\MATLAB\imdd_simulation'));
-else
-    % Linux path on the cluster
-    addpath(genpath('/work_beegfs/sutef391/imdd_simulation'));
-end
-
-% Quiet the ambiguous CET warning (best is to set TZ in sbatch; see below)
-warning('off','MATLAB:datetime:AmbiguousTimeZone');
-
-% How many workers?
-cpus = str2double(getenv('SLURM_CPUS_PER_TASK'));
-if ~isfinite(cpus) || cpus < 1, cpus = max(1, feature('numcores')); end
-
-% Use a per-job, node-local JobStorageLocation to avoid stale locks on $HOME
-% Prefer $TMPDIR if your cluster provides it, else tempdir().
-tmpbase = getenv('TMPDIR');
-if isempty(tmpbase), tmpbase = tempdir; end
-jsl = fullfile(tmpbase, sprintf('matlab_jobstorage_%s_%s', ...
-    getenv('USER'), getenv('SLURM_JOB_ID')));
-if ~exist(jsl,'dir'); mkdir(jsl); end
-
-% Configure the local cluster explicitly and start the pool
-c = parcluster('local');
-c.NumWorkers = cpus;
-c.JobStorageLocation = jsl;
-
-p = gcp('nocreate');
-if isempty(p) || p.NumWorkers ~= cpus
-    if ~isempty(p), delete(p); end
-    p = parpool(c, cpus);           % avoids the “queued” state
-end
-fprintf('parpool up with %d workers; JobStorage=%s\n', p.NumWorkers, c.JobStorageLocation);
-
-
-
-
-% result filename (timestamp + optional job id)
-t      = datetime('now','TimeZone','local','Format','yyyyMMdd_HHmmss');
-jobid  = getenv('SLURM_JOB_ID'); if isempty(jobid), jobid = 'nojid'; end
-host   = getenv('HOSTNAME'); if isempty(host), host = 'localhost'; end
-
-% Output directory depends on platform
-if ispc
-    output_root = fullfile('C:\Users\Silas\Documents\MATLAB\Datensätze\FWM_2025\');
-else
-    output_root = '/work_beegfs/sutef391/results_WDM';
-end
-if ~exist(output_root,'dir'), mkdir(output_root); end
-
-% Build filename
-t     = datetime('now','TimeZone','local','Format','yyyyMMdd_HHmmss');
-jobid = getenv('SLURM_JOB_ID'); if isempty(jobid), jobid = 'nojid'; end
-host  = getenv('HOSTNAME');     if isempty(host),  host  = 'localhost'; end
-
-fname = sprintf('WDM_%s_%s_%s.mat', char(t), host, jobid);

projects/WDM/getDispersionVector.m

@@ -0,0 +1,32 @@
+function dispersion_vector = getDispersionVector(N, D, ref_zdw, randomize_ZDW, randomkey)
+% s.MATLAB version of the Python generator shown above.
+% Returns an N×1 vector (ps/(nm·km)).
+%
+% D is the nominal dispersion magnitude. For D>0 the link is segmented with
+% alternating sign (+D, -D, +D, …). For D==0 it is flat (0) except for
+% ZDW randomization. The ZDW detuning is ~N(0, 2 nm) around 1310 nm and is
+% converted to dispersion via 0.09 ps/(nm·km) per nm.
+
+    % constants (matching the Python code)
+    meanLambda_nm = 1310;        % center wavelength
+    sigma_nm      = 2;           % ZDW sigma
+    Dslope        = 0.07;        % ps/(nm·km) per nm detuning
+
+    % random ZDW-induced dispersion offset
+    if randomize_ZDW
+        rng(randomkey, 'twister');
+        rand_zdws_nm = meanLambda_nm + sigma_nm .* randn(N,1);
+        rand_D = (rand_zdws_nm - ref_zdw) .* Dslope;   % ps/(nm·km)
+    else
+        rand_D = zeros(N,1);
+    end
+
+    % nominal segmented pattern (match Python intent; keep length N)
+    if D > 0
+        base = (-1) .^ ((0:N-1).');   % +1,-1,+1,-1,...
+    else % D == 0 (or anything else)
+        base = ones(N,1);
+    end
+
+    dispersion_vector = base .* D + rand_D;  % ps/(nm·km)
+end