;this script will generate the xraylib file containing the compton profiles as well as the second derivatives for future interpolation PRO comptonprofiles COMMON xraylib full_elec_conf = [2,$ 2,2,4,$ 2,2,4,4,6,$ 2,2,4,4,6,6,8,$ 2,2,4,4,6,6,8,$ 2,2,4,4,6,$ 2,2,4,$ 2$ ] ;read in correct biggs electronic configurations biggs_elec_config = [] OPENR,lun,'biggs_electron_configuration.dat',/get_lun line='' WHILE(NOT EOF(lun)) DO BEGIN READF,lun,line line=STRTRIM(line,1) IF (strmid(line,0,1) NE '#') THEN BEGIN biggs_elec_config = [biggs_elec_config,PTR_NEW(DOUBLE(STRSPLIT(line,/EXTRACT)))] ENDIF ENDWHILE FREE_LUN,lun FOR i=0,N_ELEMENTS(biggs_elec_config)-1 DO BEGIN IF (TOTAL(*(biggs_elec_config[i])) NE i+1) THEN BEGIN PRINT,'Error in biggs_electron_configuration.dat for element: ',i+1 STOP ENDIF ENDFOR data_template={element:0, N:0, UOCCUP: PTR_NEW(), UBIND: PTR_NEW(),pz: PTR_NEW(), total:ptr_new(), partial:ptr_new(),total_secderiv: PTR_NEW(),partial_secderiv: PTR_NEW() } nstructs=0L line='' data_reader=0 OPENR,lun,'ComptonProfiles.dat',/GET_LUN WHILE(NOT EOF(lun)) DO BEGIN READF,lun,line line=STRTRIM(line,1) SWITCH STRMID(line,0,1) OF '#': BEGIN CASE STRMID(line,0,3) OF '#S ': BEGIN IF (data_reader EQ 1) THEN BEGIN ; PRINT,'nstructs: ',nstructs ; PRINT,'N_ELEMENTS(my_data)',N_ELEMENTS(my_data) ; PRINT,'my_data[nstructs].N-2',my_data[nstructs].N-2 ; PRINT,'N_ELEMENTS(*(my_data[nstructs].partial))/(my_data[nstructs].N-2)',N_ELEMENTS(*(my_data[nstructs].partial))/(my_data[nstructs].N-2) *(my_data[nstructs].partial) = ALOG(REFORM(*(my_data[nstructs].partial), my_data[nstructs].N-2,N_ELEMENTS(*(my_data[nstructs].partial))/(my_data[nstructs].N-2),/overwrite)) *(my_data[nstructs].total) = ALOG(*(my_data[nstructs].total)) *(my_data[nstructs].pz)=ALOG(*(my_data[nstructs].pz)+1) ;calculate second derivatives my_data[nstructs].total_secderiv = PTR_NEW(DERIV(*(my_data[nstructs].pz),DERIV(*(my_data[nstructs].pz),*(my_data[nstructs].total)))) my_data[nstructs].partial_secderiv = PTR_NEW(DBLARR(my_data[nstructs].N-2,N_ELEMENTS(*(my_data[nstructs].partial))/(my_data[nstructs].N-2))) FOR i=0,my_data[nstructs].N-3 DO BEGIN (*(my_data[nstructs].partial_secderiv))[i,*] = DERIV(*(my_data[nstructs].pz),DERIV(*(my_data[nstructs].pz),(*(my_data[nstructs].partial))[i,*])) ENDFOR ;adapt partial comptonprofiles print,'n elec_config:',N_ELEMENTS(*(biggs_elec_config[nstructs])) print,'elec config:',*(biggs_elec_config[nstructs]) print,'n uoccup:',N_ELEMENTS(*(my_data[nstructs].UOCCUP)) print,'uoccup:',*(my_data[nstructs].UOCCUP) IF (N_ELEMENTS(*(biggs_elec_config[nstructs])) NE N_ELEMENTS(*(my_data[nstructs].UOCCUP))) THEN BEGIN j=0 partial_new=DBLARR(N_ELEMENTS(*(biggs_elec_config[nstructs])),N_ELEMENTS(*(my_data[nstructs].pz))) partial_secderiv_new=DBLARR(N_ELEMENTS(*(biggs_elec_config[nstructs])),N_ELEMENTS(*(my_data[nstructs].pz))) for i=0,N_ELEMENTS(*(my_data[nstructs].UOCCUP))-1 DO BEGIN print,'i:',i print,'j:',j print,'elec config j:',(*(biggs_elec_config[nstructs]))[j] print,'uoccup',(*(my_data[nstructs].UOCCUP))[i] my_zeroes = DBLARR(1,N_ELEMENTS(*(my_data[nstructs].pz))) WHILE ((*(biggs_elec_config[nstructs]))[j] EQ 0.0) DO BEGIN partial_new[j,*] = my_zeroes j = j+1 ENDWHILE IF ((*(biggs_elec_config[nstructs]))[j] LT (*(my_data[nstructs].UOCCUP))[i]) THEN BEGIN partial_new[j,*] = (*(my_data[nstructs].partial))[i,*] partial_new[j+1,*] = (*(my_data[nstructs].partial))[i,*] partial_secderiv_new[j,*] = (*(my_data[nstructs].partial_secderiv))[i,*] partial_secderiv_new[j+1,*] = (*(my_data[nstructs].partial_secderiv))[i,*] j = j+2 ENDIF ELSE BEGIN partial_new[j,*] = (*(my_data[nstructs].partial))[i,*] partial_secderiv_new[j,*] = (*(my_data[nstructs].partial_secderiv))[i,*] j = j+1 ENDELSE ENDFOR ptr_free,my_data[nstructs].partial ptr_free,my_data[nstructs].partial_secderiv ptr_free,my_data[nstructs].UOCCUP my_data[nstructs].partial = PTR_NEW(TEMPORARY(partial_new)) my_data[nstructs].partial_secderiv = PTR_NEW(TEMPORARY(partial_secderiv_new)) my_data[nstructs].N=N_ELEMENTS(*(biggs_elec_config[nstructs]))+2 my_data[nstructs].uoccup = biggs_elec_config[nstructs] ENDIF nstructs++ data_reader=0 ENDIF IF (nstructs EQ 0) THEN BEGIN my_data = [REPLICATE(data_template,1)] ENDIF ELSE BEGIN my_data = [my_data,REPLICATE(data_template,1)] ENDELSE splitted=STRSPLIT(line,/EXTRACT) my_data[nstructs].element=FIX(splitted[1]) print,'Element: ',my_data[nstructs].element END '#N ':BEGIN splitted=STRSPLIT(line,/EXTRACT) my_data[nstructs].N=FIX(splitted[1]) END '#UO':BEGIN splitted=STRSPLIT(line,/EXTRACT) my_data[nstructs].UOCCUP=PTR_NEW(FIX(splitted[1:my_data[nstructs].N-2])) ;analyze electronic configuration ;construct kissel array ;elec_config=[] ;for i=K_SHELL,Q3_SHELL do begin ; temp_value= ElectronConfig(my_data[nstructs].element,i) ; elec_config = [elec_config, temp_value] ;endfor ;nonzero = WHERE(elec_config GT 0.0) ;elec_config = elec_config[0:nonzero[-1]] ;IF (N_ELEMENTS(elec_config) EQ N_ELEMENTS(*(my_data[nstructs].UOCCUP))) THEN BEGIN ; PRINT,'Electronic configuration OK for '+STRING(my_data[nstructs].element) ;ENDIF END '#UB':BEGIN splitted=STRSPLIT(line,/EXTRACT) my_data[nstructs].UBIND=PTR_NEW(DOUBLE(splitted[1:my_data[nstructs].N-2])) END '#L ': BEGIN data_reader=1 END ELSE: ENDCASE BREAK END '0': '1': '2': '3': '4': '5': '6': '7': '8': '9': BEGIN IF (NOT data_reader) THEN BEGIN PRINT,'data_reader error: aborting' PRINT,'Last string was: ',line return ENDIF values=DOUBLE(STRSPLIT(line,/EXTRACT)) IF (NOT PTR_VALID(my_data[nstructs].pz)) THEN BEGIN my_data[nstructs].pz = PTR_NEW([DOUBLE(values[0])]) my_data[nstructs].total = PTR_NEW([DOUBLE(values[1])]) my_data[nstructs].partial = PTR_NEW([DOUBLE(values[2:my_data[nstructs].N-1])]) ENDIF ELSE BEGIN *(my_data[nstructs].pz) = [*(my_data[nstructs].pz),DOUBLE(values[0])] *(my_data[nstructs].total) = [*(my_data[nstructs].total),DOUBLE(values[1])] *(my_data[nstructs].partial) = [*(my_data[nstructs].partial),DOUBLE(values[2:my_data[nstructs].N-1])] ENDELSE END ENDSWITCH ENDWHILE *(my_data[nstructs].partial) = ALOG(REFORM(*(my_data[nstructs].partial), my_data[nstructs].N-2,N_ELEMENTS(*(my_data[nstructs].partial))/(my_data[nstructs].N-2),/overwrite)) *(my_data[nstructs].total) = ALOG(*(my_data[nstructs].total)) *(my_data[nstructs].pz)=ALOG(*(my_data[nstructs].pz)+1) ;calculate second derivatives my_data[nstructs].total_secderiv = PTR_NEW(DERIV(*(my_data[nstructs].pz),DERIV(*(my_data[nstructs].pz),*(my_data[nstructs].total)))) my_data[nstructs].partial_secderiv = PTR_NEW(DBLARR(my_data[nstructs].N-2,N_ELEMENTS(*(my_data[nstructs].partial))/(my_data[nstructs].N-2))) FOR i=0,my_data[nstructs].N-3 DO BEGIN (*(my_data[nstructs].partial_secderiv))[i,*] = DERIV(*(my_data[nstructs].pz),DERIV(*(my_data[nstructs].pz),(*(my_data[nstructs].partial))[i,*])) ENDFOR IF (N_ELEMENTS(*(biggs_elec_config[nstructs])) NE N_ELEMENTS(*(my_data[nstructs].UOCCUP))) THEN BEGIN j=0 partial_new=DBLARR(N_ELEMENTS(*(biggs_elec_config[nstructs])),N_ELEMENTS(*(my_data[nstructs].pz))) partial_secderiv_new=DBLARR(N_ELEMENTS(*(biggs_elec_config[nstructs])),N_ELEMENTS(*(my_data[nstructs].pz))) for i=0,N_ELEMENTS(*(my_data[nstructs].UOCCUP))-1 DO BEGIN my_zeroes = DBLARR(1,N_ELEMENTS(*(my_data[nstructs].pz))) WHILE ((*(biggs_elec_config[nstructs]))[j] EQ 0.0) DO BEGIN partial_new[j,*] = my_zeroes j = j+1 ENDWHILE IF ((*(biggs_elec_config[nstructs]))[j] LT (*(my_data[nstructs].UOCCUP))[i]) THEN BEGIN partial_new[j,*] = (*(my_data[nstructs].partial))[i,*] partial_new[j+1,*] = (*(my_data[nstructs].partial))[i,*] partial_secderiv_new[j,*] = (*(my_data[nstructs].partial_secderiv))[i,*] partial_secderiv_new[j+1,*] = (*(my_data[nstructs].partial_secderiv))[i,*] j = j+2 ENDIF ELSE BEGIN partial_new[j,*] = (*(my_data[nstructs].partial))[i,*] partial_secderiv_new[j,*] = (*(my_data[nstructs].partial_secderiv))[i,*] j = j+1 ENDELSE ENDFOR ptr_free,my_data[nstructs].partial ptr_free,my_data[nstructs].partial_secderiv my_data[nstructs].partial = PTR_NEW(TEMPORARY(partial_new)) my_data[nstructs].partial_secderiv = PTR_NEW(TEMPORARY(partial_secderiv_new)) my_data[nstructs].uoccup = biggs_elec_config[nstructs] my_data[nstructs].N=N_ELEMENTS(*(biggs_elec_config[nstructs]))+2 ENDIF nstructs++ FREE_LUN,lun ;reading in seems to work ;so does calculating the second derivatives ;write them to a file OPENW,lun,'../comptonprofiles.dat',/GET_LUN FOR i=0,nstructs-1 DO BEGIN PRINTF,lun,my_data[i].N-2,N_ELEMENTS(*(my_data[i].pz)) PRINTF,lun,FIX(*(my_data[i].UOCCUP)) PRINTF,lun,*(my_data[i].PZ) PRINTF,lun,*(my_data[i].TOTAL) PRINTF,lun,*(my_data[i].TOTAL_SECDERIV) FOR j=0,my_data[i].N-3 DO BEGIN IF ((*(my_data[i].UOCCUP))[j] LE 0.0) THEN CONTINUE ;FOR k=0,N_ELEMENTS(*(my_data[i].pz))-1 DO BEGIN PRINTF,lun,(*(my_data[i].PARTIAL))[j,*] ;ENDFOR ENDFOR FOR j=0,my_data[i].N-3 DO BEGIN IF ((*(my_data[i].UOCCUP))[j] LE 0.0) THEN CONTINUE ;FOR k=0,N_ELEMENTS(*(my_data[i].pz))-1 DO BEGIN PRINTF,lun,(*(my_data[i].PARTIAL_SECDERIV))[j,*] ;ENDFOR ENDFOR ENDFOR FREE_LUN,lun END